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ABSTRACT 


The tactical implications of submarine acoustic radiation and UNDEX-survivability 
have motivated the development of an advanced machinery cradle which will provide shock 
and vibration isolation of the submarine internals, thereby minimizing the resulting acoustic 
radiation. The cradle space frame must be designed and optimized for both minimum 
shock/vibration bi-directional transmissibility and minimum total cradle weight. Frequency 
domain structural synthesis (structural modification and substructure coupling), is applied 
to the cradle design. The method addresses static and complex dynamic problems in 
structural design analysis, and allows the direct analytic treatment of specialized equipment, 
such as frequency-dependent visco-elastic isolators. 
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1. INTRODl'CTION 


The design of complex structural dynamic systems requires the building of detailed 
mathematical models with which to predict static and dynamic response. Most commonly, 
the finite element (FE) method is used to generate structure system matrices with which 
dynamic response can be calculated. While the FE method currently provides the best 
means of predicting response for complex strucmral systems, the time required to assemble 
the system matrices and to process them for the calculation of dynamic response can be 
prohibitive. Therefore, the use of the FE method for performing design analyses often 
precludes the performance of numerous design analyses in the search for an optimal 
design. This is especially true when a FE based analysis is to be used in conjunction with 
advanced design techniques such as optimization. The iterative process of modeling the 
system and analyzing the model to determine the system performance is the design-analysis 
cycle. 

The traditional design-analysis cycle consists of the following process. A designer 
builds a FE model which best represents the system. The system model is the complete 
system structure, for example, a submarine hull and an internal machinery support cradle is 
modeled as one structure. The definition of the FE model yields system matrices, which 
include stiffness, mass, and less commonly damping. The numerical generation of the 
system matnces is referred to as the assembly phase. At this point, loads are apphed and 
responses, static and/or dynamic, are calculated. The calculation of system response is 
referred to as the solution phase. The responses are then used to calculate stresses and 
strains in the model. These calculations are referred to as the post-processing phase. The 
solution phase is the most costly in terms of time and computing resources. 




Based on the acceptability of the displacements, stresses, and strains calculated, the 
design or system model may have to be changed in the interest of improving the response 
characteristics of the initial or follow on design. For example, a high stress which is 
unacceptable may exist at a certain location in the design. The designer decides that if a 
particular alteration is made to the design, the stress response will become within tolerable 
levels. Traditionally, this alteration requires a repeat of the assembly, solution, and post¬ 
processing phase of the analysis, a cost and time intensive procedure which limits the 
number of design re-analyses that can be accomplished. Since the re-analysis is time 
consuming, the optimal design is abandoned for a final design which is less than optimal. 

Therefore, with the intent of accelerating the design process and lowering the attendant 
costs, new methodologies for assembling and modifying system models is presented. The 
new method replaces all three of the FE analysis phases with a single computationally 
efficient calculation. The method to be described herein, generally referred to as frequency 
domain structural synthesis [Refs. 1,2,3], is directed specificaUy at drastically reducing the 
time required to perform a design analysis cycle. This capability for rapid re-analysis makes 
structural synthesis ideal for use in advanced automated design environments, such as in 
conjunction with optimization codes. 



II. FINITE ELEMENT EORMl LAI ION 


The finite element method used for comparison with the solution obtained from 
frequency domain structural synthesis is based on Lagrange's equation of motion [Ref. 4|. 
Various types of elements are used in modeling of structural systems, including for 
example, plate, shell, and beam elements. Our discussion will be limited to beam elements 
expenencing combined bending and axial deformations. We are using beam elements to 
demonstrate the methodology because the beam element allows fora manageable system of 
equations and matrices that are easily handled by a personal computer. The theory remains 
valid for all types of elements and is unaffected by the complexity of shell or plate 
elements. Beam elements that are subject to bending and axial deformation have three joint 
displacements at each end of the beam element The beam element has six generalized 
coordinates and six degrees of freedom (DOF), which yields a mass, stiffness, and 
damping matrix for the beam element of size (6 x 6). The beam element is shown in 
Figure 1. 



Figure 1. Beam Element with Coordinate and Nodal Orientation 

Each node has a set of coordinates, axial (x), lateral (y), and rotational (0), however 
not limited to three tX)F. Elements can be modeled with six IX)F per node. 
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elements are 




The derivation presented in reference(41 assumes that the axial forces associated with the 
axial joint displacements (x) have only a negligible effect on the shape functions associated 
with the joint displacements (y) and (0). With this assumption the mass and stiffness 
matrices are derived. 

The elemental stiffness and mass matrices are 



where the terms E is Young’s Modulus, I is the area moment of inertia, 1 is the elemental 
beam length in inches, y is the weight density, and r is the radius of gyration. The 
elemental matrices are partitioned in the following way: 

^ x e, „ 02 


The damping matrix is usually impossible to determine analytically and is typically 
determined experimentally. Here damping is applied to the system model by one of three 
ways. The three methods generally used are: 

Type (1): Proportional structural damping of the form: 


[C] = a[A:] + p[A/l 


(1) 



Type (2): Proportional viscous damping of the form: 

(C| = af/fl ,2, 

Type i3l; Frequency-dependent viscous damping of the form: 

[C’l = [C,e-'“] (3) 

Type (1) damping is used in adding damping to structural elements, and Types (2) and (3) 
are used in adding damping to vibration isolators which are a combination of springs and 
dampers; adding proportional damping to just the isolators constitutes non-proportional 
damping for the whole structure. 

The equation of motion of these finite elements can be written in terms of their joint 
displacements as 

K.|«,L+[4i«.LMHkl = UL (4) 

where («J ^ = axial, lateral, or angular joint displacements 

[m]^ - mass matrix of element 
[c], = damping matrix of element 
[/t], = stiffness matrix of element 
\f\^= joint forces and moments 

Since the elements vary in orientation with respect to the system axis, the elemental 
mass and stiffness matrices must be transformed into global coordinates. The 
transformation of the damping matrix is neglected, since the elemental damping matnx is 
modeled as a function of the transformed elemental mass and stiffness matrix. A method 
for relating local joint displacements of each element to the global system displacements 
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must be incorporated. This method is referred to as a coordinate transformation. The 
elemental mass and stiffness matrices after transformation are in the following form 
[Ref. 4|. 


' 1400 + I56S- 
-16CS 
y/ ‘ -22/S 

420' 70C- + 54S- 
; I6CS 

13/S 


-I6CS 

140S- + 156C’ 
22/C 
16CS 

70S- + 54C- 
-I3ZC 


-22/S 70C- +54S- 

22/C 16CS 

4/- -13/S 

-13/S 140C-+156S^ 
13/C -16CS 

-31- 22/S 


16CS 

70S- + 54C- 
13/C 
-16CS 

140S- + 156C- 
-22IC 


- CS-12CS -6/S - -IC--12S- 


CS+12CS -6/S 


CS-12CS S^ + 12C^ 

-6/S 




iiY 


CS+12CS 6/S 


-f- CS + 12CS -f-ls--12C- 6/C 

I'-J 

^ 6/s -bic 21 - \ 

f-1 C- + 12S^ f-1 CS-12CS 6/S 


CS + 12CS 

-bis 


- 12C^ -6/C 


CS-12CS 


flY.: 


K 12C- 


-6/C 


Noting that: 

1 = elemental beam length in inches 
E = Young's modules in psi 
I = Area Moment of Inertia in 
r = radius of gyration in inches 
Y = mass density per unit length in Ih ■/ /«’ 
C = cos a 


S = sin a 



Once the Finite elements are transformed to global coordinates, the elements are 
assembled to generate the global mass, stiffness, and damping matrices. The equation of 
motion for the modeled system in global coordinates is 

|M||i/|>[Cll«M^l|«i = (5) 

where \u\ = axial, lateral, or angular global joint displacements 
\M\ = global mass matnx 
[C|= global damping matrix 
[^]= global stiflfiess matrix 
I F] = global joint forces and moments 

The following example will demonstrate how the global mass and stiffness matrices are 
generated. Consider the structural system modeled with two beam elements and having no 
boundary conditions shown in Figure 2. 

1,2,3 ^ 4,5,6 / 7,8,9 

: ~ I I 

1 2 3 

Figure 2. A Beam Modeled Using Two Elements 

El yl 

For this example y =1.0 Ibs/in, ^^ = 1.0 lb s , I = 1.0 in and r=1.0 in.. Since the 

beam elements lie horiTontally along the x axis, the angle a = 0. The global mass and 
stiffness matrices are generated from the assembly of the elemental matrices. The elemental 


matrices: 





■ 1 0 0-1 

0 12 6 0 



0 - 12-6 0 
0 6 2 0 



M40 0 0 70 0 

0 156 22 0 54 

0 22 4 0 13 

70 0 0 140 0 

0 54 13 0 156 

_ 0 -13 -3 0 -22 


0 

-13 

-3 

0 


Referring to Figure 2. the lower bold type numerals represent the node numbering and the 
upper numbers represent the beam nodal coordinates. Coordinates 1, 2. 3. 4. 5, 6, 7. 8, 
and 9 are respectively -r|.y..0..;c„y,,9,.x,,y,, and 9,. Remembering how the matrices are 
partitioned and noting that coordinates 4, 5 and 6 of beam element 1 are the same 
coordinates of beam element 2 and thus are shared. The two elemental matrices are 
assembled together through the shared coordinates. Figure 3 shows the beam element 
arrangements. 


1.2.3 


1 



^ 8.9 


3 


Figure 3. Beam Elements with Node and Global Nodal Coordinate Numbering 


For discussion purposes only, the stiffness matrices will be demonstrated, since the 
mass and damping matrices are generated in the same manner. The elemental stiffness 
matrices are in the following form. 




8j 


2 3 


5 6 


5 6 7 




The two matnces are combined by adding shared nodal coordinates, this process is 
detemiined by the element connectivity. Figure 3 shows that element 1 is coupled to 
element 2 through global nodal coordinates 4. 5 and 6. These global coordinates are the 
combination of local ctxirdinates x.. \\. 0, of element 1 and .t.. v,. 0, of element 2. The 
resulting matrix is a (9x9) global matrix represented by [All, The size of the global rnatnx 
IS the number of nodes times the DOF. The global stiffness rnatnx (Mj is shown below 
with the numbers installed, take special note to the shared coordinates which are additive. 


■ 1 

0 

0 

-1 

0 0 

0 

0 0 ■ 

0 

12 

6 

0 

-12 6 

0 

0 0 1 

0 

6 

4 

0 

-6 2 

0 

0 0 

-1 

0 

0 

2 

0 0 

-J 

0 0 

\K\= 0 

-12 

-6 

0 

24 0 

0 

-12 6 , 

0 

6 

2 

0 

0 8 

0 

-6 2 

0 

0 

0 

-1 

0 0 

1 

0 0 

0 

0 

0 

0 

-12 -6 

0 

12 -6 

L 0 

0 

0 

0 

6 2 

0 

-6 4 ( 


The shared coordinates are demonstrated by looking at the 3x3 partition, rows 4 through 6 
and columns 4 through 6. After the global matrix is generated, the boundary conditions are 
applied. Boundary conditions are determined by coordinate restraints. If a coordinate is 
restrained then the row and column corresponding to that coordinate are deleted. For 
example, if in Figure 2, the left end had been fixed, displacements for global coordinates 1, 
2. and the slope of coordinate 3 are zero, and therefore the rows and columns 
corresponding to these coordinates would be deleted resulting in a (6x6) [ K | matrix. 

The derivation of the equation for a second order linear structural system described in 


the frequency domain is presented below. The differential equation of motion for a second 
order linear structural system is written as 





+ CX + kx = f sinflA. 


The solution to equation (6), which is the total system response is 

= iX",. + iX", (7) 

where is the real or homogeneous solution and is the particular solution. We 
consider only steady state harmonic excitation, therefore the particular solution is used as 
the total solution. The solution is assumed to have the form 

X = X^=Xe^\ ( 8 ) 

Taking the first and second derivatives of X and substituting into equation (6) yields 

{-kfmX + jClcX + kxy^^ = . (9) 

Dividing both sides by and rearranging equation (9) gives the equation of motion as 

{k-Cl-m+jQc)x = F. (10) 

Writing equation (10) in matrix form gives the equation for second order linear 
structural systems described in the frequency domain. 


[[x]-finwhyfafi)i]uiHFi 



where the vector |.tl is the set of generalized responses in the global coordinate system. 
The vector (Ft contains generalized global forces and moments. [F| and \M\ are 
symmetnc. real valued and of order n. The damping matrix [O Q)] is in general, frequency 
dependent, but here is modeled as a linear proportional combination of the mass and 
stiffness matrices. Equation (6) is generally written in compact form as 

[Z(0'i]|.r| = (Ft tl 2 ) 

where the matrix [Z(n)] is called the system impedance matrix. Equadon (12) is the system 
impedance reladonship and represents the dynamic response of the system. The impedance 
matrix is the dynamic stiffness of the system. The static case is when 0 = 0 and then the 
system impedance matrix is just the stiffness matrix [ F ]. 

The impedance matrix for the assembled beam (Figure 2) is 

[Z\ = [K]-D-[M]. (13) 

Equation (13) is reduced because damping [CtO)] is neglected in this example. [Z(0)] is 
calculated over the frequency band of interest, where O is the frequency band of mterest in 
rad/sec. 

The frequency response function (FRF) matrix for the assembled beam is 


[//(0)] = [Z(0)f‘. 


(14) 





The FRF matrix allows the calculation of the steady state harmonic response amplitude j Xl 
resulting from a harmonic force amplitude {Fj. The frequency response relation is 
determined by matrix inversion of equation (12), which yields 

ui =[//(a)]|Fl. (151 

Any element of the frequency response matrix is defined as the dynamic response of 
motion coordinate i due to a unit harmonic generalized force acting on motion coordinate j. 
The FRF matrix can be used to represent information about displacements, velocities, 
accelerations, stress, or strains. For example, if a structure is excited at nodal coordinate 5, 
then //,, is the complex amplitude of the response at nodal coordinate 1 due to a unit 
harmonic excitation at nodal coordinate 5 at some frequency Q, of interest. 

A typical frequency response function plot is shown in Figure 4. The peaks shown in 
Figure 4 occur at the frequency of peak response. The relationship between the natural 
undamped frequency, the damped natural frequency, and the frequency of peak response is 
shown on the following page. The plot shows at what frequencies the structure will have 
maximum responses and enables the designer to redesign the structure so that the system 
will have small responses in the frequency bandwidth of interest. 
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Frequency Hz 


Figure 4. Typical Frequency Response Function Plot 


Generally there are three distinct frequencies of interest, the undamped natural 
frequency, the damped natural frequency, and the frequency of peak response. These 
frequencies are related by the modal damping factor. The amplitude of the response of a 
forced vibration can become very large when the frequency of the excitation approaches 
one of the natural frequencies of the system. This condition where the excitation frequency 
IS the same as one of the natural frequencies is referred to as resonance. When a system 
vibrates at resonance, the attendant stresses and strains have the potential of causing 
structural failure. A structural system will have a maximum response when the frequency 
of excitation is near the undamped natural frequency. If the system has no damping, then 
the maximum response will occur at the undamped natural frequency. The undamped 
natural frequency is a function of the system mass and stiffness and is analytically 
expressed as the solution to the eigensystem 


= |01 


(16) 






Every real structural system has an infinite number of natural frequencies and mode shapes. 
The finite modeling of the structural system yields a finite number of eigensolutions or 
mode shapes and eigenvalues or natural frequencies depending on how many degrees of 
freedom the structural system is modeled with. Each eigenvector has a corresponding 
eigenvalue or natural frequency. However, all systems inherently have some degree of 
damping and the relationship that relates the damped natural frequency to the undamped 
natural frequency is 


to^ = io„.yi-C,' (17) 

where C, is the modal damping factor for mode i. The damped natural frequency is slightly 
lower than the undamped natural frequency and a typical damping factor for structural 
systems is 0.2. The frequency of peak response is the frequency of excitation where the 
response of the system is maximum. The analytical relation that relates the frequency of 
peak response to the undamped natural frequency is determined by taking the derivative 
with respect to j equation (18) and setting it equal to zero. 



and substituting back into equation (18) yields 


1^1 1 

J[i-tb^]-’ + [2(;,wf ■ 

Now performing = 0 

uix: 

- :^|(l - (b - )■’ + (2C,tb )■' |'^[ -4u)(l - (I) -) + 8C,-(1)] = 0 


(19; 


(20) 


and knowing for equation (20) to equal zero, the numerator must equal zero. Setting the 
numerator equal to zero 


2tb(l -u)’)-4C,'tb = 0 


and simplifying 


(b 


^ =1-2C/ 


( 22 ) 


and solving for to , which is the frequency of peak response yields 


(23) 


It is important when designing a structural system that the excitation frequency is not close 
to these frequencies or failure of the structural system is likely to occur. 




III. FREQUENCY DOMAIN STRUCTLTRAE SYNTHESES 


The theory presented herein is taken directly and exclusively from references [1,2. and 
3|. The purpose of this thesis is to explore the application of this previously developed 
theory to the analysis of a submarine cradle structure. 

Frequency domain structural synthesis was first presented in 1939 and has evolved to 
the latest formulation, which was published in Journal of Sound and Vibration (1991) 
[Ref. Ij. The most recent formulation of the theory is a new method for analyzing 
structural systems [Ref. 3]. This method handles all types of structural models and is more 
efficient and cost effective compared with traditional finite element solution procedures. 
Frequency domain structural synthesis refers to substructure coupling and structural 
modification using frequency response function data. The previously developed 
formulation for structural synthesis, Ref. [3], is applicable to the static and dynamic 
structural analysis of direct coupling of substructures, indirect coupling of substructures, 
modification of substructures, and constraint application. The theory allows the synthesis 
of displacements, velocities, accelerations, stresses, and strains. 

An important feature of the frequency domain formulation is the arbitrary and exact 
model order reduction possible when performing a synthesis. A Finite element method 
(FEM) when applied to practical problems typically generates between KP to 10^ degrees- 
of-freedom (DOF). The frequency domain formulation allows, as a minimum, only those 
DOF of interest to be included in the analysis. This feature is in fact the reason for the high 
computational efficiency of the method. Using one of the numerical examples presented in 
the section “Numerical Examples, ” the computing time required for a frequency domain 
synthesis can be compared with the same analysis using traditional finite element (EE) 
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procedure. Referring to Example (6) the following count of floating point operations 
(FLOPS) shows the efficiency of the frequency domain method: 

FEM direct assembly: Time - 25876 sec or 431.3 mins 
FI .OPS - 1.49 X 10^ 

FRF synthesis: Time - 1167 sec or 19.45 mins 
FLOPS-517.2 X I06 

This clearly demonstrates that synthesis by FRF is more efficient and better suited for the 
re-analysis of complex structures with large numbers of DOF. Moreover, the savings in 
time grows with increasing model size. 

There are two major classifications of structural synthesis. These classifications are 
coupling and modification and each classification can be viewed as direct or indirect. 
Coupling is defined as the joining of two separate substructures to form one structure and 
modification is defined as the creation of a new load path in an existing structure. An 
example of coupling is the coupling of a submarine hull and the machinery support cradle. 
The hull is modeled as one substructure and the cradle is modeled as another substructure. 
We want to join these two substructures together to create one complete structural system. 
This process is known as structural coupling. As an example of the use of structural 
modification, an analysis of the complete cradle structural system shows that a certain 
element has unacceptable stresses. By installing an additional support, the stresses become 
acceptable. This process of changing the structure is known as structural modification. 
Indirect coupling is the joining of structures with the introduction of an intermediate or 
interconnecting structural element, referred to as an interconnection impedance element or 
impedance patch [Ref. 1], Direct modification can be viewed as the application of a 
constraint equation to a given structural model; direct coupling is simple substructure 
synthesis [Ref. 1], The theory is unique in that it allows any linear structural element to be 
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used as an interconnection impedance, for example a spnng and viscous damper may be 
installed between two elements of a structure. The synthesis is performed at each frequency 
of interest, which makes possible the efficient treatment of frequency dependent properties, 
like the properties in the spring and viscous damper. Frequency domain structural synthesis 
allows changes to a finite element model without reassembly of the mass, stiffness, or 
damping matrices. 

A. GENERALIZED FREQUENCY RESPONSE 

The derivations presented here are taken exclusively from References 1, 2. and 3, The 
derivations are reproduced with more intermediate steps leading to the final operative 
equations. We begin the development with the previously derived formulation for a second 
order linear structural system described in the frequency domain. 

The differential equation of motion for a second order linear structural system is 
written as 


mx + cx + kx = F sin Clt. (6) 

The solution to equation (6), which is the total system response is 

X = (7) 

where X^ is the real or homogeneous solution and X^ is the particular solution. We 
consider only steady state harmonic excitation, therefore the particular solution is used as 
the total solution. The solution is assumed to have the form 

X = X^=X€^. ( 8 ) 
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Taking the first and second derivatives of X and substituting into equation (6) yields 

(-n-' mX + jClrX + kX)e^'‘‘ = Fe'"". (9) 

Dividing both sides by and rearranging equation (9) gives the equation of motion as 

{k-Q-m-<-jac]X= F. nO) 

The dynamic stiffness of the structural system is known as the system impedance which is 
wntten as 

Z{Ci)^k-Q.-m + jnc. (24) 

The static system stiffness k is determined by the case where D = 0 and the system 
impedance is 

Z(Q)=k. (25) 

The matrix notation for the structural system impedance is 

[Z(n)]|xl = i/1. (26) 

The system impedance matrix [Z(D)] is both complex valued and frequency dependent. 
The general equation for the frequency response structural model is found by taking the 
matrix inversion of equation (26) and is indicated as 
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=[«(a)]{/fn)|. 


[ 21 \ 


|x} and 1/1 are vectors of complex valued generalized response and excitation coordinates 
at a specific frequency fl. and IW] is the frequency response function (FRF) matrix 
evaluated at the frequency fi. In general an element of the FRF matrix is defined by talcing 
the partial derivative of |a-| with respect to |/|. Referring to equation (24) and writing the 
equation for jr, we get 


-q = //,,/ + A + W,,/, + ■■■ + ( 28 ) 

and taking the partial derivative of equation (28), the general form for an element of the 
FRF matrix is 


H. = 


dx^ 

% 


(29) 


and is defined as the partial derivative of the ith generalized response coordinate with 
respect to the jth generalized excitation coordinate. 

There are other types of frequency response which are classified by the type of 
coordinates involved. For example, strain-force and stress-force frequency response are 
defined as 


|E|=1«']I/1 |al=Kll/|. (30,31) 

The difference between equation (27) and the two equations (30,31) is the FRF matrix [ H\ 
contains displacement-force information in equation (27) and strain/stress-force information 
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in equations (30.31). The general element of the strain and stress FRF is determined in the 
same manner as the displacement-force FRF. The general elements are defined as 




where e, and o, are complex valued strains and stresses at coordinates i at a specific 
frequency D 

Flere we will show the development of the frequency response function in the modal 
coordinate system. We start with the differential equation of motion for a second order 
linear structural system in physical coordinates. 


where we assume F(r) is of the form and the vector jf) is the set of force 

amplitudes. Now we apply the linear transformation 

= (35) 

to equation (34) where the vector j.r| is the set of physical coordinates to be transformed, 
the vector |^) is the set of modal coordinates, and the matrix [<1>] is the set of mass 
normalized normal mode shapes. Pre multiplying the transformed equation by and 
using the relation [<t>] = [/] yields 

[/]|<7l+[ 2C,w, jl^lJ w: ||.7l=[4>f{Fl ={3^1- (36) 
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Rewriting equation (36), the differential equation of motion in the modal coordinate system 
is 


+ 2^,10;^ + = .|(r). 


(37) 


Since we assumed steady state harmonic forces, equation (34), the modal forces are also of 
the form |J) = and the solution is assumed as a steady state harmonic modal 

response of the form {<?( = |cV \e^'. Taking the first and second derivative of {< 7 ! and 
substituting back into equaf '7) and simplifying yields 


to; -Tl-[/l+d 


2C,w, 




(38) 


Rewriting equation (38) 


to: -n'-rya 2 C,io, 


bl = {j| 


(39) 


and solving equation (39) in terms of the modal response yields 



(o;-n^ryn2C,to, 


Pl- 


(40) 
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To transform equation (40) back to the physical coordinate system, we use the 
transformation of modal force. [<t>] jf'l = |Jl. and the transformation of modal 
coordinates. |-r( = [cl>]{,V ), to substitute back into equation (40) and simplify. The 
resulting equation in physical coordinates is 


\x] = m 


d!,- - a' + 7x12;;,to, 



Remembering the general form of the frequency response, equation (27). the frequency 
response function [//(fl)] in terms of the system modal information is 


[Mn)] = [ci>] 


to ■ - n- + 7x12;;,to, 




(42) 


and any specific element of H is given by 


Hi, = 


I 


to; -X4' + 7X12C,to/ 


(4.4) 


B. MATRIX PARTITIONING 

First we will define the classification of coordinates. Figure 5 represents two 
substructures A and B that will be joined together by merging coordinates 2 with 3 and 6 
with 7. These coordinates are referred to as connection coordinates and are denoted by the 
subscript "c". By the definition just stated, the connection coordinates for substructure A is 
2 and 6, hkewise the connection coordinates for substructure B are 3 and 7. Internal 
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coordinates, denoted by the subscript "i," are all the remaining coordinates not directly 
involved in the substructure coupling. In Figure 5, the internal coordinates for substructure 
A are 1 and 5 and the internal coordinates for substructure B are 4 and 8. The set of all the 
physical coordinates are denoted as coordinate set "e". If one structure is involved, then the 
coordinate set "e" contains only the connection and internal coordinates for that structure. If 
two or more substructures are involved, then the coordinate set "e" contains all the 
coordinates for all the substructures. The mathematical representation is e = i U c. 

nr:n 

5 6 7 8 

A B 

Figure 5. Structural Model with Internal and Connection Coordinates 

Referring to the general equation for frequency response, equation (27), and writing it 
in matrix form with coordinate partitioning as 

Hi; am 

where x) and fi are a set of generalized responses and excitations at the internal coordinates 
and xc and fc are a set of generalized response and excitation at the connection coordinates. 
One of the special features about the frequency response is that in addition to response 
information, we can also determine other information at the same time, for example, 
stresses. We can append a set of stress coordinates and then equation (44) becomes 
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The stress coordinates will allow the direct calculation of synthesized system stress. 

The generalized excitations are partitioned into internal and external excitations. 
Referring to Figure 5, and looking at the internal coordinates, for example coordinate 1, it 
is obvious that the only force possible on this coordinate is an externally applied force. 
Since the internal coordinates do not participate in synthesis, there are no coupling forces 
present on internal coordinates. The connection coordinates, for example, coordinate 2 in 
Figure 5, may experience both externally applied forces and coupling forces which are 
established through synthesis. Therefore 

( 46 ) 

and by definition of the internal coordinates 

/=/“ (47) 


Introducing equations (46 and 47) into equation (45) allows for the expansion of equation 
(45) as 



25 







where the asterisk superscript denotes a synthesized quantity due to the fact that we have 
introduced the forces of synthesis, Note that with the introduction of equations (46 
and 47). a redundant equation, the fourth row of equation (48), has been appended. Using 
the definition of the set "e", equation (48) is written in the new condensed form 



(49) 


where \f^ j = [/'" ] . The vector fe is externally applied forces which may exist at 

all physical coordinates, and the vector fc is the coupling forces present only at the 
connection coordinates. 


C. STRUCTURAL MODIFICATION ANT) INDIRECT SUBSTRUCTURE 

COUPLING 

In this section we will develop the governing equation for structural modification and 
indirect substructure coupling. As previously defined, structural modification is the creation 
of redundant load paths within a structure, and indirect coupling is the creation of new load 
path with a structural element between uncoupled substructures. Indirect coupling and 
structural modification are confined to connection coordinates. There is only one restriction 
enforced for these processes. The structural change used for modification or 
interconnection impedance used for indirect coupling must be described by the following 
equation 


) = -[ K{a} - M{Q) + ;c(n )]{x;) (SO) 

or 
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where the negative sign shows that the reaction is on the structure to be modified or 
substructures to be coupled. Equation (50) defines the transformation of forces which is 
used to transform equation (49), The transformation which operates on equation (49) is 


oiM-l 

1/.J Lo -zJU;j 


Substituting equation (52) into equation (49) yields 


(53) 


then performing matrix multiplication and simplifying equation (53), the resulting form is 
given by the relationship 






-£z| 




Extracting the third row of equation (54) 


(54) 


(55) 


and rearranging equation (55) 
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[/^W, Z]|.r;) = [//,Jl/,j ,56) 

and solving equation (56) in terms of j.r’} yields 

|.r;)=[/.W..,zr'[//JU). (57) 

Extracting the second row of equation (54) 

and substituting equation (57) into equation (58) 

= [W«]|/,l - l«.c][Z]([/ ^ (59) 

now, using the known relation that 

k‘i=[w:]ui (60) 

and substituting equation (60) into equation (59) gives the following relation 

[w„]‘ui=[w„]ui -\hmi +«„zr'[//„]ui ( 61 ) 

and rearranging equation (61) and setting it equal to zero 

l//«]‘l/e} Iz][/ + //„zr'[//„]lx) = (01 (62) 
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and rewnting equation (62) 


1[««1'" «-,.ZriW.,.|]l/J = |0[. (63) 

Since by definition {/. j ’= |0|, then 

[[«,.r -[ W.„ 1 + [ W.. I Z\[l + = [0] (64) 

and solving equation (64) in terms of [ W„] yields 

= [^..1 - (65) 

Now we will simplify the third term of equation (65). Extracting the following portion 

factoring the inverse term 

[Z][(Z-‘r//„)z]‘' 

and applying the identity (a<>) = a''h^' 

[Z]\Z'{Z-'^ Hj'\ 
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and then simplifying yields 


Substituting the above simplified portion back into equation (65) and performing the same 
process on the first row of equation (54) yields the final operative equation for structural 
modification and indirect coupling 



Note for the static case when = 0. the impedance matrix [Z] = [AT] and [Af] is a singular 
matrix which is not invertible, therefore equation (66) is not valid and a form of the 
equation which does not require the matrix inversion of [Z] must be used. The following 
equation is for the static case when = 0. 



Terms on the right side of equations (66 and 67) are pre-synthesized values and the left 
hand side is the synthesized values. The matrix [Z] describes the modification to be made 
for structural modification and can be negative if the modification to be made is the removal 
of a structural modification, or it describes the new load path between two structures for 
indirect coupling. The quantity [//„J allows for the direct calculation of stress due to 
externally applied loads in the synthesized structure. Stress frequency response could be 
replaced with strain or other structural frequency response. 
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I). SI BSTRLCTURE COI PI.ING AND CONSTRAINT IMPOSITION 

In this section, the development of the theory of direct substructure coupling using 
boolean mapping matrices is shown. A formal discussion of the mapping matnx is 
presented in the next section. The development of this theory also applies to constraint 
imposition. Substructure coupling involves the joining of two or more separate 
substructures where constraint imposition involves one structure, the coupled structure, 
Constraint imposition is the application of two conditions on the synthesized connection 
coordinates. The first condition being force equilibrium where the summation of forces on 
a coordinate are equal to zero and the second condition being compatibility where the 
displacement of the synthesized coordinates are equal to zero. Compatibility is interpreted 
as the connection coordinate from the first substructure most have the same displacement as 
the connection coordinate from the second substructure in order for them to be merged as a 
single coordinate. 

We extract the third row from equation (49) which is shown here again for reader 
convenience. 


[ ^ 1 








r' 

= 

[\ 

(49) 

kl 





The third row of equation (49) is 

(68) 

We construct the conditions for equilibrium and compatibility to be imposed on the 
connection coordinates. Figure 6 shows two connection coordinates from two 
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substructures and Figure 7 shows the equilibrium and compatibility conditions applied to 
the merged connection coordinate of the synthesized structure. 


c c 

• ◄-► • 

Figure 6. Connection Coordinates from Two Substructures 





Figure 7. Merged Connection Coordinate 


Referring to Figure 7, we write the equilibrium equation for the pair of connection 
coordinates shown in Figure 6 as 


=0 (69) 

where the superscript denotes the substructure and the compatibility equation for the 
merged connection coordinates is 


;t;-x;=0. (70) 

Converting equations (69 and 70) into the general equations which will encompass all the 
connection coordinates. The general form uses the mapping matrix to relate each pair of 
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coordinates. Noting that /' = -f' and .r‘ = .r. we can write the general equations. The 
general equation for the force equilibrium is 

(71) 

where the vector |/. j represents the arbitrarily selected independent subset of the 
connection coordinates. Noting that the mapping matrix [W] must remain constant for the 
constraint imposition to hold, the general form of the compatibility equation is 

rt,.i=[Mrixj = iot (72) 

where the vector (x.j represents the compatibility for the pairs of connection coordinates. 
This vector is the zero vector. 

The transformation equations that operate on equation (49) are derived from the general 
equilibrium and compatibility equations and are of the form 



Substituting equation (73 and 74) into equation (49) 


KT [/ 0T//„ Hi! 0¥/l 

U,.J Lo //„lo mJUJ 
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noting that we are using the displacement frequency response only for derivation purposes, 
and then performing matrix multiplication and simplifying equation (75). The resulting 
form is given by equation (76). 


HJ4 1[/;l 


(76) 


Extracting the second row of equation (76) 

= (77) 

and rewriting the second terra of the left hand side 


and enforcing compatibility between pairs of connection coordinates. 


{i.l = (0| 


(72) 


equation (70) becomes 


loi = ["cX]Ul+[4.]U|- 


(78) 


Rearranging equation (78) 


[wjj/4 =-[ h , xu \ 


(79) 
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and pre multiplying both sides of equation (79) by [ W.,.] yields 


i/h-iw-xi" r'l/.i- ( 80 ) 

Extracting the first row of equation (76) 

=[««l|/ehKW]|X| ( 81 ) 

and substituting equation (80) into equation (81) 

= \h,m\ - \hMhx\kX m ( 82 ) 

and substituting the known relation 

l'r;i =[<]{/,) ( 60 ) 

for the right hand side of equation (82) 

[««rui = \HM\ U1 (83) 


and dividing both sides by (^) yields the operative equation for dtiect substructure 
coupling 

[hJ = [//„] - [//,<■ ][M][w„]'‘[ A/''J//„] (84) 
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where the terms on the right hand side are frequency response values calculated from the 
uncoupled structures and the left hand side is the frequency response values for the coupled 
system. 


Performing the same derivation presented above on the first row of equation (49), 
yields the operative equation for direct substructure coupling with coupled system stress 
response. 




E. DIRECTED GRAPHS AND MAPPING MATRICES 

The theory of direct substructure coupling requires mapping matrices to invoke the 
constraints of equilibrium and compatibility. The theory developed in the preceding section 
demonstrated how the mapping matrix represented the conditions of equibbrium and 
compatibiUty on the synthesized connection coordinates. The mapping matrix can be 
constructed from a graph which represents the connectivity that is estabbshed when 
substructures are coupled through synthesis. The general formulation of the mapping 
matrices using directed graphs presented below is taken directly from reference [3]. 


The use of equation (85) to perform substructure coupbng requires the construction 
of the mapping matrices, [M], As was developed in the pret^ng section, each column 
of [M] represents a statement of the equilibrium and compatibility which is enforced for 
each pair of connection coordinates bang coupled. We will now demonstrate that fMl 
can be constructed from a graph which is drawn to represent the connectivity to be 
estabbshed through the synthesis. 
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Figure 8. Substructure Couplings and I>irected Graphs 


Consider the coupling depicted on the left in Figure 8. Substructure "A" is being 
coupled to substructure “B,” through, say, a single pair of connection coordinates, x 
and X. Tfie coupling of this pair of coordinates creates load path “I.” To construct the 
mapping matrix for this connection “I,” we arbitrarily assign a value of ’1” to the 
connection coordinate of substructure “A” and a value 1" to the connection coordinate 
of substructure ■‘B.” The mapping matrix for this connection is 


[M] 



A" 

B" 


Considering now the more complicated coupling on the right of Figure 8, and also 
acknowledging that in general two substructures are coupled using more than one pair 
of connection coordinates, we may construct the mapping matrix. Here, the 
connections “I”, “J”, and “K” consist of more than one pam of connection coordinates 
each; these are, in general, sets of connection coordinate pairs. The mapping matrix is 


"A" 

"B" 

"C" 


where each column contains plus/minus identity matrices whose elements correspond to 
the coupling to be established between each pair of connection coordinates. For 
example, in column 2 of the above mapping matrix, all connection coordinates 
associated with substructure “A” are assign^ a “1” (i.e. HI) and they are to be coupled 
to their counterparts in substructure “C” which have been assigned a “-1” (i.e. 

The coupling of these coordinates constitutes the set of load paths denoted as “J”. 

The directed graphs and their boolean mapping matrices provide a means of 
organizing complex couplings, and also provide a framework for the computational 
implementation of the synthesis, i.e. equation (85). Of course, care must be exercised 
to insure that all matrices in equation (85) are appropriately partitioned. 


37 





An example of using directed graphs to generate the mapping matrix is presented here. 
Figure 9 shows substructure coupling and the associated directed graph for the load paths 
created when connection coordinates are coupled through synthesis 



Figure 9. Substructure Coupling Using Directed Graphs 

The upper portion of Figure 9 shows two nodes from two substructures that are to be 
coupled and the lower portion of Figure 9 shows that each node has three degrees of 
freedom which correlates to three coordinates. Each coordinate of node 1 is synthesized to 
its corresponding coordinate of node 2. The synthesis of the these coordinates creates load 
paths "A", "B", and "C”. Invoking the constraints of equilibrium and compatibility, the 
mapping matrix is constructed. Using the equilibrium equation presented earlier 

U} = lM](2| (71) 

where the vector j^) is the complete set of connection coordinates from both substructures 
and the vector j is the arbitrary selected subset of connection coordinates to be retained 

pertaining to the selected substructure. From the equilibrium equation, we get the 
relationship between the two subsets of connection coordinates, = -f^. If we arbitrarily 
select the connection coordinates from substructure 1 as our set to retain, then we assign a 
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1 to those coordinates and from the relationship shown above, we assign a -1 to the 
connection coordinates of substructure 2. The mapping matrix tor the system in Figure 9 is 
determined using equation (71) 


'A- "B" 'C- 
■ 1 0 0" 


hi 


( 86 ) 


where 

■ I 0 0] 

0 1 0 


0 -1 0 I 
1,0 0 -IJ 

The upper partition of equation (86) corresponds to arbitrarily selected coordinates of 
substructure 1 and the lower partition corresponds to connection coordinates of 
substructure 2. The mapping matrix relates how theses coordinates are connected. 

F. MODIFICATION AND INDIRECT COUPLING USING MAPPING 
MATRICES 

This section will show the development of the operative equation of synthesis for 
indirect substructure coupling and structural modification. There are two classes of 
synthesis for which mapping matrices are used The first class is direct substructure 
coupling which was discussed in Section C and the second class is for indirect substructure 
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coupling and structural modification. This class of synthesis again uses interconnection 
impedance to synthesize two substructures or modify an existing structure. The mapping 
matrix contains the connectivity information corresponding to the equilibrium of the 
interconnection impedance and the equilibrium of the modification. The interconnection 
impedance for this method of synthesis has the requirement that the structural element used 
as an interconnection impedance must be described without mass terms. The 
interconnection impedance is a function of stiffness and damping. This method is well 
suited for the synthesis of visco-elastic isolators between substructures. 

Visco-elastic isolators are modeled as a combination of a spring and dash pot damper. 
The isolators are treated as having proportional viscous dampers or as having frequency 
dependent viscous dampers. A special note here is that adding proportional damping to just 
the isolators constitutes non-proportional damping for the complete synthesized strucmre. 

We begin the derivation with the description of the structural system, equation (49). 


Ct 1’ 


hJ 


-cj = 



.^J 





The transformation matrices which operate on equation (49) and lead to the operative 
equation for indirect coupling and modification using mapping matrices are 



and 


[ctT f/ 0 Olfaj’ 

X, [ = 0 / 0 Li 

LJ LO 0 


( 88 ) 
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The impedance introduced in equation (87) is a reduced system impedance that is massless 
and is of the form 


|z|=[Mf[Z]([,Wl")\ 


(89) 


Using the displacement frequency response of equation (49) and substituting the 
transformation equations, equation (87 and 88) into equation (49) 


[■/ 0 T W . Hi 1 

0 1|/1 

LO W„lo 

1 -MzJliJ 


and simplifying equation (90) yields 


J/1 

U,J 


(91) 


Extracting the second row of equation (91) 


{l\ = (92) 


and rewriting equation (92) yields 

+ |(W"//„A/z]{jf,l' = 1 • (93) 


Equation (93) is rewritten so that the left hand side is a product of sums 




[/-/w^//,,wz||:r,j'=[w^w,„]jx). (941 

Pre multiplying both sides of equation (94) by [/+ ' and simplifying yields 

|x,|’ = [/ - . (95) 

Extracting the first row of equation (91) 

(96) 

and substituting equation (95) into equation (96) 

M' =K11/,1 + K1[A/]|z]/+ (97) 

and using the definition of the frequency response 

\<\-K\U\ (60) 

to substitute into the left hand side of equation (97) yields 

= [Hj - [//„][A/][z]/ + (98) 

Noting that 


42 



we will simplify the third term of equation (94). Extracting the following portion 






factoring the inverse term 


|zIlr‘-H„)z|' 


and then simplifying yields 


[i-’ * H,,]"- 


Substituting the above simplified portion back into equation (98) yields the final operative 
equation, equation (99), for indirect synthesis and modification using mapping matrices 



(99) 


where the terms on the right hand side are frequency response values calculated from the 
uncoupled structures and the left hand side is the frequency response values for the coupled 
system. 

Performing the same derivation presented above on the first row of equation (49), 
yields the operative equation for indirect substructure coupling and modification using 
mapping matrices with coupled system stress response. 



(100) 
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I\. NUMERICAL KXAMPI.KS 


The following numerical examples are provided to give a detailed explanation for each 
type of synthesis. The results of each example are presented graphically and are compared 
with the traditional finite element method (FEM) solution. 

Three types of damping that are addressed in the numerical examples are: 

Type (1): Proportional structural damping of the form: 

[C] = a[/fl + p[A/l (1) 

Type (2): Proportional viscous damping of the form: 

[C] = a[/!:] (2) 

Type (3): Frequency-dependent viscous damping of the form: 

[C] = [c„e-‘'“] (3) 

Type (1) damping is used in adding damping to a substructure, and Types (2) and (3) are 
used in adding damping to the isolators which are a combination of spring and dampers. 
The system impedance matrix Z for these damping types are: 

Type (1): [Z(n)] = [/f] - + iC], where [C] = a[/f] + p[W]. (100) 

Type (2): [Z(n)] = [AT] - n^M] + ^Cl, where [C] = af/C]. (101) 

Type (3): [Z(n)] = [/s:]-n-[A/] + /l(an)], where [C] - [Qe-”"]. (102) 
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A. EXAMPLE (1): DYNAMIC INDIRECT COCPEING 


Consider the structures shown in the following figures. The structure shown in Figure 
1.1 will be directly assembled by the finite element method in order to compare a traditional 
calculation of the frequency response with that synthesized from the substructures shown 
in Figure 1.2. 


Figure 1.1. Structure Analyzed Using Traditional FE Procedures. 



Figure 1.2. Synthesis of Structure 


The total structure shown in Figure 1.1 is synthesized from Structure 1 and Structure 2 
through the interconnection impedance Z or “new load path.” For this example, the 
following beam parameters will be used: 

Young's Modulus E= 30.0 » 10* psi 
Area moment of inertia I = 0.1666 * 10'^ in'" 

Cross-sectional area A = 0.2 in^ 

Weight density WTD = 0.2832 Ibf/in^ 

2 percent proportional structural damping a = 0.02 
Beam element lengths = 24 in 

Proportional structural damping is applied only to structures 1 and 2, and the 
interconnection impedance Z is undamped. The damping applied in this example was 
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arbitrarily selected. The synthesis method is not limited to proportional structural damping, 
any arbitrary linear frequency dependent dampmg can be used. Referring to Figure 1.2, the 
system of structures 1 and 2 and the interconnection impedance z are synthesized in the 
frequency domain to yield exact results as the FEM direct assembly method. The general 
synthesis equation for dynamic indirect coupling is 

- [Hecp ' * ( 66 ) 

Note again that structures 1 and 2 have proportional structural damping and the 
interconnection impedance Z is undamped. The general procedure for performing the 
synthesis is as follows. 

The mass and stiffness matrices [/k] and [ Af] for the three substructures (including the 
middle beam analytically treated as an interconnection impedance) are generated using 
traditional FEM. Since each structure is comprised of only one beam element, the elemental 
matrices with boundary conditions applied are the substructure global matrices. The 
impedance matrix is calculated for each structure as 

. (103) 

Note that [/f,] and [A^] are complex-valued and [M,] is real-valued. The FRF matrix [H\ 
for structures 1 and 2 is calculated by inverting the impedance matrix [Z\. We now have 
[Wj], [//,] and [zj. Referring to the general synthesis equation provided above, the 
matrices [W„], [W„], and [//„] ate generated by assembling [//,] and [W,] by 

appropriate partitioning. [//„] is the combination of [//,] and [W,] and is partitioned by 
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internal and connection coordinates. Referring to Figure 1.2. after the boundary conditions 
are applied, coordinates 4, 5. and 6 are renumbered 1. 2, and 3 respectively for structure 1. 
Structure 2 is unaffected since the boundary conditions remove coordinates 4. 5, and 6 and 
coordinates 1. 2. and 3 remain the same. The impedance z is unchanged. [//.J is 
partitioned in the following manner 




HiU) H(U) 
Mc,i) H(c,cl 


where the subscript “i” represents the set of internal coordinates and "c" represents the set 
of connection coordinates. A more detailed representation is 




h 

c, 


h 


0 


0 


0 


0 



0 


0 

c. 

0 


0 

H,{c„cA 


In this representation “il" denotes the internal coordinates of structure 1 and “cl” denotes 
the connection coordinates of structure 1. The same principle follows for “i2” and “c2” 
relating to structure 2. The partitioning for Hec, Hcc, and Hce are 








0 




q[«,(c„q) 0 

H^{c2,c 
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^ cT//,(c.O 0 |//,(c-,x) 0 

c,[ 0 0 

In this example there are no internal coordinates. The connection coordinates for structure 1 
are (1, 2, 3) and for structure 2 are (1. 2, 3). With the appropriate partitioning complete, 
the synthesis can now be performed. Using the indirect coupling relation 

[«.,]. ( 66 ) 

structure 1 is synthesized to structure 2 through the "new load path" Z. [f/,J is the 
synthesized FRF relation representing the exact dynamics of the total structure. Finally the 
FRF relation is calculated over the frequency range 0.1 - 65 Hz and plotted in the Figures 
which follow. 



Figure 1.3. Plot of Hee (2,2) from Synthesis 







Frequency Hz 

Figure 1.4. Plot H (2,2) from Traditional FE Calculation 

Figures 1.3 is a plot of the synthesized [//„| matrix, element (2,2), and Figure 1.4 is 
the same FRF element calculated using the traditional FE procedure. The FRF element 
plotted in both figures corresponds to coordinate 5 of Figure 1.1, a lateral motion 
coordinate. Notice both plots are identical, demonstrating that the synthesis procedure 
provides an exact solution for the synthesized system dynamics. The figures show the first 
four damped natural frequencies. The FRF plots show the magnitude of the response at 
coordinate 5 due to a unit excitation of varying frequency at coordinate 5. 

B. EXAMPLE (2): DYNAMIC DIRECT COUPLING 

Consider the following figures. The structure shown in Figure 2.1 will be directly 
assembled using FEM for the purpose of comparing with the results obtained by 
synthesizing structures 1 and 2 of Figure 2.2. 
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Figure 



2.1. Hull-Cradle Structure Analyzed by Traditional FE Techniques 



Structure 2 will be coupled to structure 1 at coordinates 10, 11, 12, 13, 14, 15, 31. 32, 33, 
34, 35, and 36. These coordinates are the connection coordinates and the remaining 
coordinates are internal coordinates. Coordinates 1, 2, 3, 13, 14, 15, 16, 17, 18 19, 20, 
and 21 of structure 2 are connection coordinates and the remaining are internal. The 
following beam element data will be used: 

Young's Modulus E = 30.0 « 10‘ psi 
Area moment of inertia 1 = 0.02083 in" 

Cross-sectional area A = 1 in 
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Weight density WTD = 0.2832 ibf/in' 

Proportional structural damping (1 %) a = 0.01 
The proportional structural damping was arbitrarily selected and is applied to both 
structures. The general equation for dynamic direct coupling is 

( 84 ) 

In this equation, M is the boolean mapping matrix which is used to establish the 
connectivity between the two substructures for synthesis. The mapping matrix is 
determined by the connectivity i.e. what is connected to what and by imposing the 
equilibrium and compatibility relations associated with each pair of coordinates. We can 
define the mapping matrix by j/) = [Wl|/,.j. Where j/,) is a vector of aU the comiection 

coordinates of both structures and j/,. | is the arbitrarily selected independent subset of the 

connection coordinates relating to one of the substructures. We have selected the 
connection coordinates of structure 1 as the arbitrary subset of connection coordinates. The 
mapping matrix [A/] is a matrix of size (24 x 12) and is depicted as; 
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We will calculate the FRF matrix [H] tor both substructures. First the [^] and [M\ 
matrices are generated for each substructure. [AfJ and [A",] are both complex since 
proportional structural damping was applied to both structures; again this damping is 
arbitrary, [/f. ] and [/fj are of the form [a:] = [A: + jaK]. We next form the impedance 
matrix for each substructure. The impedance matrix is of the form [Z] = [A"] - Q’[A/]. With 
the impedance matrix generated for each substructure, the FRF matrix FI can be calculated 
by inverting the impedance matrix. This process is done at each frequency of interest. 
These FRF matrices are required in order to couple the two structures together to form the 
structure in figure 2.1. Referring to the synthesis equation above, the matrices [//„]. 
[W„.], [W,], and [W.,] are formed by combining [//,] and [W.) by appropriate 
partitioning. The partitioning is shown below. 
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Referring to Figure 2.2, "i I” denotes the set of internal coordinates of structure 1 which are 
1. 2, 3, 4. 5, 6. 7. 8, 9, 16. 17, 18. 19. 20, 21. 22, 23. 24. 25, 26. 27. 28, 29. and .30. 
"cl” denotes the set of connection coordinates of structure 1 which are 10, 11, 12, 13, 14, 
15. 31, 32. 33, 34, 35. and 36. “i2" denotes the set of internal coordinates of structure 2 
which are 4. 5. 6, 7, 8, 9. 10, 11. and 12, "c2’' denotes the set of connection coordinates 
of structure 2 which are 1. 2. 3. 13, 14. 15. 16. 17, 18. 19. 20. and 21. With the 
appropriate partitioning complete, the synthesis of structure 1 to structure 2 can be 
performed using the direct coupling relation 

[«.,]. (84) 

[//„] is the synthesized FRF relation which is the combination of both structures. The 
synthesis is done over the frequency range of interest and plotted in Figure 2.3. The 
frequency range for this example was 0.1 to 10.0 Hz. Figure 2.4 is the solution from 
traditional FE calculations included for direct comparison. Both plots are identical. 
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Figure 2.4. Plot of H (8,8) from Traditional FE Calculations. 

Figures 2.3 and 2.4 are the plots of the FRF at element (8,8) from the synthesized and FE 
[//] matrices. This element corresponds to the lateral motion coordinate 8 of Figure 2.1. 
Notice both plots are identical and both show the first seven damped natural frequencies. 
The plots show the magnitude of the response of unit amplitude at coordinate 8 due to a 
unit excitation at varying frequency at coordinate 8. As the frequency of excitation 
approaches the damped natural frequency, the response approaches infinity. 

C. EXAMPLE (3): STOUCTLIRAL MODIFICATION (RFJMOVAL OF A 
BEAM ELEMENT) 

Consider the following figures. Figure 3.1 depicts a combined hull-cradle structure 
which will be directly assembled by traditional FE procedures. Note that the structure in 
Figure 3.1 has asymmetric reinforcing trusses. The synthesis methodology will be used to 
arrive at the structural configuration shown on the left of Figure 3.1 by removing the beam 
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shown on the right of Figure 3.2. The FRF calculated from the FE model (Figure 3.1) will 
be compared with that calculated using synthesis. 



Figure 3.1. Final Hull-Cradle Configuration 



Structire 1 Structure 2 

Figure 3.2. Synthesis Used to Remove a Beam Element 


Referring to Figure 3.2, structure 1 will be modified by removing the beam, structure 2. 
located between nodal coordinates 10, 11, 12, 43, 44, and 45. The following beam element 
data will be used; 

Young's Modulus E = 30.0 >< 10* psi 
Area moment of inertia I = 0.02083 in'* 

Cross-sectional area A = 1 in 
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Weight density WTD = 0.2832 Ibf/in' 

1 percent proportional structural damping a = 0.01 
The proportional structural damping was arbitrarily selected and is applied to both 
structures. The general equation for dynamic indirect coupling/modificabon is 

= [«..] - [«..][«■, - (66) 


Note that the sign in the term -Z"'] is opposite from that in the original indirect 
coupling equation. This is because we are removing the beam element from the structure 
instead of synthesizing it to the structure. The first step is to generate the [/f] and [M] 
matrices for structure 1 and structure 2. The [/f| matrices for both structures are complex 
since proportional damping was applied. They are of the form [ Al] = [Af + jaX\. Next we 
form the impedance matrices for each structure, [Z] = [Af] - This method requires 

the calculation of the FRF matrix [ W] only for the structure to be modified, structure 1 of 
Figure 3.2, The impedance and the FRF matrices are calculated at the frequency of interest. 
Once the and impedance matrices are generated, we are ready to partition the FRF 
matrix. The matrices [w„], [W ,.], and [w ,] are formed by partitioning [ wj . The 
partitioning is shown below. 
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The connection coordinates for structure 1 are 10. 11. 12. 43. 44. and 45. The rest are all 
treated as internal coordinates. With the appropriate partitioning of [//,] completed, the 
removal of the beam from the structure can now be completed by using the correct form of 
the indirect coupling relation mentioned above. [//^J is the synthesized FRF relation 
which reflects the removal of structure 2 from of structure 1. This modification is 
calculated over the frequency range of interest and plotted in Figure 3.3. The frequency 
range for this example was 0.1 to 7.0 Hz. Figure 3.4 is the solution from the traditional 
FE procedure and is provided to allow direct comparison of the two solutions. Both plots 
are identical. 



Figure 3.3. Plot of Hee(l 1,11) as Calculated Using Synthesis 
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Figure 3.4. Plot of H( 14.14) Calculated Using Traditional FE Procedures 

Figures 3.3 and 3.4 are the plots of the FRF corresponding to the lateral motion coordinate 
14 of Figure 3.1. A special note here is that the element (14,14) of the FRF generated by 
FEM is the coordinate 14, which corresponds to the element (11,11) of the FRF generated 
by the indirect coupling relation. The reason for this is because of the partitioning. [ H ,,) is 
partitioned with internal coordinates first followed by the connection coordinates. Care is 
required here to ensure the coordinate of interest is actually being used. Notice both plots 
are identical and show the first six damped natural frequencies. The plots show the 
magnitude of the response at coordinate 14 due to a unit excitation at varying frequency at 
coordinate 14. As the frequency of excitation approaches the damped natural frequency, the 
response approaches infinity. 
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D. EXAMFLK (4): STRUCTURAL MODIFICATION (ADl)THON OF A 
BEAM) 

Consider the following figures. The FRF for the structure shown in Figure 4.1 will be 
calculated by traditional FE procedures to compare with that calculated using the synthesis 
procedure to add the beam element, as shown in Figure 4.2. 



Figure 4.1. Hull-Cradle Structure Analyzed by Traditional FETechmques. 



Strurture 1 Struaure 2 

Figure 4.2. Synthesis is Used to Add the Beam Element 
Referring to Figure 4.2, structure 1 wiU be modified by adding the beam, structure 2, at the 
nodal coordinates 10, 11, 12, 43, 44. and 45. The following beam element data was used: 
Young's Modulus E = 30.0 x lO" psi 
Area moment of inertia 1 = 0.02083 in'* 
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Cross-sectional area A = 1 in 
Weight density WTD = 0.2832 Ibf/in' 

Proportional stnjctural damping (1%) a = 0.01 
The proportional structural damping was arbitrarily selected and is applied to both 
structures. The general equation for dynamic indirect coupling/modification is 

^ ]" J ■ (66) 


The first step is to generate the [Af] and [A/] matrices for structure 1 and structure 2. The 
[AT] matrices for both structures are complex since proportional damping was applied. 
They are of the form |Ar] = [Al+ jaK\. Next, impedance matrices are formed for each 
structure as [Z] = [^] - n'[M]. This method requires the calculation of the FRF matrix 
[H\ only for the structure to be modified, structure 1 of Figure 4.2. The impedance and the 
FRF matrices are calculated at the frequency of interest. Once the FRF and impedance 
matrices are generated, partition of the FRF matrix is required. The matrices [ , [ W,, ]. 

and [//,,] are formed by partitioning [//,]. The partitioning is shown below. 
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The connection coordinates for structure 1 are 10, 11, 12, 43, 44, and 45. The rest are all 
treated as internal coordinates. With the appropriate partitioning of [//,] completed, the 
synthesis of the beam to the structure can now be completed by using the correct fon: of 
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the indirect coupling relation mentioned above. [W,„] is the modified FRF relation which 
is the combination of structure 1 and the added element, structure 2. The synthesis is 
performed over the frequency range of interest and plotted in Figure 4.3. The frequency 
range for this example is 0.1 to 8.5 Hz. Figure 4.4 is the solution from a traditional FE 
calculation for direct comparison of the two solutions. Both plots are identical. 



Figure 4.3. Plot of Synthesized FRF Element Hee{8.8) 
m 200 - - - 

I 100- 
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8 -100 , - - - 
“ -200-- - 

^0123456789 
Frequency Hz 

Figure 4.4. Plot of H (8,8) Calculated Using Traditional FE Procedures. 
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Figures 4.3 and 4.4 are the plots of the FRF corresponding to the lateral motion coordinate 
8 of Figure 4,1. A special note here is that the element (8.8) of the FRF generated by FEM 
corresponds to the coordinate 8, as does the element (8.8) of the FRF generated by the 
indirect coupling relation. This is different from the previous example. The reason for this 
is because of the partitioning. \H^^\ is partitioned with internal coordinates first followed 
by the connection coordinates. Care is required here to ensure the coordinate of interest is 
actually being used. Notice both plots are identical and show the First six damped natural 
frequencies. The plots show the magnitude of the response at coordinate 8 due to a unit 
excitation at varying frequency at coordinate 8. As the frequency of excitation approaches 
the damped natural frequency, the response approaches infinity. 

E. EXAMPLE (5); INDIRECT COUPLING WITH ISOLATORS 

Consider the following figures. The FRF for the structure shown in Figure 5.1 wiU be 
calculated using traditional FE procedures to compare with the FRF calculated using the 
synthesis method. The synthesis will combine the various components shown in Figure 
5.2. In this figure, the hull model (structure 1) will be coupled to the cradle model 
(structure 2). Note that this example demonstrates that the synthesis procedure easily and 
exactly treats problems with non-proportional damping, a truly unique feature of the 
methodology. 
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Figure 5.1. Traditional FE Procedures are Used to Calculate FRF for the Combined 
Hull-Isolator-Cradle Structural System. 



Figure 5.2. Total Hull-Isolator-Cradle System is Synthesized from Components. 


Referring to Figure 5.2, structure 1, structure 2, and four spring-damper isolator sets will 
be synthesized together to form the system in Figure 5.1. Each isolator set consists of three 
spring-damper isolators, one for each connection coordinate. The connection coordinates 
for structure 1 are 10, 11, 12, 13, 14, 15, 31, 32, 33, 34, 35, and 36. The remaining 
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coordinates are internal coordinates. Coordinates 1. 2, 3. 13, 14, 15, 16. 17. 18 19. 20. 
and 21 of structure 2 are connection coordinates and the remaining are internal. For this 
structural synthesis method, the spring-damper isolators are treated as a lumped system 
(with no physical dimensions) installed at the connection coordinates. The connection 
coordinates do not merge into one but are joined by way of the isolators. The following 
beam element data will be used: 

Young's Modulus E = 30.0 » 10^ psi 
Area moment of inertia 1 = 0.02083 in'* 

Cross-sectional area A = 1 in 
Weight density WTD = 0.2832 Ibf/in' 

Proportional structural damping (2%) a = 0.02 
Proportional viscous damping (2%) P = 0.02 
Isolator spring constant k = 25 Ibs/in 

The proportional structural damping was arbitrarily selected and is applied to both 
structure 1 and 2 of Figure 5.2. The proportional viscous damping used for the damper in 
the isolator is arbitrary and is not limited to being proportional but could be any frequency 
dependent function. For our example the isolator is of the analytic form [ k + j Opk] where 
j = Recalling the impedance relation [Z(fl)] = [Af] - + yn[C] , [C] is the 

proportional viscous damping, [P^]. The operative equation for indirect coupling with 
mapping matrices is 


(99) 

where [z| = and [h,,] = [mY[H„\m]. 
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Note that [z] reduces to [/|(k + jOPk) and its size is ( 12 x 12 ). The boolean mapping 
matrix [Mils determined the same way as explained in example two. The connection 
coordinates for structure 1 and structure 2 are listed above. We can define the mapping 
matrix by |/.| = [A^]|/^.j. Where {/.j is a vector of all the connection coordinates of both 

structures and |/,. j is the arbitrarily selected independent subset of the connection 

coordinates relating to one of the substructures. We have selected structure 1 as the 
arbitrary subset of connection coordinates. The mapping matrix \M\ is a matrix of size (24 
x 12) and is 


The FRF matrix [//] for both substructures is required. First the [Af] and [/V/| matrices are 
generated for each substructure, [/f,] and [/f,] are both complex since proportional 
structural damping was applied to both structures, again this damping is arbitrary. [ K, ] and 
[/f,] are of the form [/(■] = [/(■ + joK]. We next form the impedance matrix for each 
substructure. The impedance matrix is of the form [Z] = [/f] - Q^[ Af]. With the impedance 
matrix generated for each substructure, the FRF matrix H can be calculated by inverting the 
impedance matrix. This process is done at each frequency of interest. Now with the FRF 
matrix for each substructure calculated, we are ready to synthesize the two structures and 
isolators together to form the structure in Figure 5.1. Referring to the synthesis equation 
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above, the matrices [/7„], [W,,,]. and [//^,] are formed by combining [//J and [//,] 

by appropriate partitioning. The partitioning is shown below. 
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Referring to Figure 5.2, “il” denotes the set of internal coordinates of structure 1 which 
include 1, 2, 3, 4, 5, 6, 7, 8, 9, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 ,29, 
and 30, “cl” denotes the set of connection coordinates of structure 1 which include 10, 11, 
12, 13, 14, 15, 31, 32, 33, 34, 35, and 36, “i2” denotes the set of internal coordinates of 
structure 2 which include 4, 5, 6, 7, 8, 9, 10, 11, and 12, “c2” denotes the set of 
connection coordinates of structure 2 which are 1, 2, 3, 13, 14, 15, 16, 17, 18, 19, 20, 
and 21. With the ap propriate partitioning complete, the synthesis of structure 1 to structure 
2 can be performed using the indirect coupling relation 
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(99i 


1^-1'=[W.J - 

Structure 1 is synthesized to structure 2 by the isolators or load paths described by [zj. 
[// ^] is the synthesized FRF relation which is the combination of both structures and 
isolators. The synthesis is done over the frequency range of interest and the response is 
plotted in Figure 5.3. The frequency range for this example was 0.1 to 8.0 Hz. Figure 5.4 
is the solution from traditional FE calculations for direct comparison. Both plots are 
identical. 



Frequency Hz 

Figure 5.3. Plot of Hee(8,8) for Synthesized System. 
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Figure 5.4. Plot of H(8.8) from Traditional FE Calculations. 


Figures 5.3 and 5.4 are the plots of the FRF at element (8,8) of the synthesized \H^^\ and 
traditional FE [//] matrices, respectively. This element (8.8) corresponds to the lateral 
motion coordinate 8 of Figure 5.1. Notice both plots are identical and show the first five 
damped natural frequencies. The plots show the magnitude of the response at coordinate 8 
due to a unit excitation at varying frequency at coordinate 8. As the frequency of excitation 
approaches the damped natural frequency, the response approaches infinity. 

F. EXAMPLE (6): INDIRECT COUPLING WITH FREQUENCY 

DEPENDENT ISOLATORS 

In this example, we demonstrate the capability of synthesizing components with 
frequency dependent properties. Specifically, we will repeat the preceding example using 
isolators that have frequency dependent damping. Consider the following figures. The FRF 
for the structure shown in Figure 6.1 wUl be calculated using traditional FE procedures to 
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compare with the FRF calculated using the synthesize method. The components to be 
synthesized are shown in Figure 6.2. 


Figure 6.1. Traditional FE Procedures are Used to Calculate FRF for the Combined 
Hull-lsolator-Cradle Structural System 




Sprng vfd Qimper Isolator 

Figure 6.2. Total Hull-Isolator-Cradle System is Synthesized from Components. 

In Example (5), the isolators were treated as having proportional viscous damping. This 
example will use viscous damping which is frequency dependent. The methodology is the 
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same as in the previous example and will not be repeated here. The discussion here will 
focus on the only difference which is the damping applied to the isolator. Referring to the 

general impedance relation [ZCQ)] = [^f] - O’[A/] + yO[C] . [C) is now a function of fi. 

The equarion is now of the form 

[z(a:i] = [A:j-n-’[/v/] + /i(c(n)]. (i 04 ) 

The damping applied to the isolator was arbitrarily selected as an exponential decay 
dependent on frequency. The form of the function used is 

C=C/^ 0) 

where C„ = A: = 25 lb • s/in and a = 0.1. 

The damping function is plotted in Figure 6.3. 



The reduced impedance in this example is now of the form [z] = [/](k + jDke'‘“) and the 
size of the matrix is (12 x 12). 
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Structure 1 is coupled to structure 2 with isolators described by [z]. [W,.^] is the 
synthesized FRF relation which is the combination of both structures and isolators. The 
synthesis is done over the frequency range of interest and plotted. Figure 6.4. The 
frequency range for this example was 0.1 to 8.0 Hz. Figure 6.5 is the solution from 
traditional FE calculations provided for direct comparison. Both plots are identical. 



Frequency Hz 

Figure 6.4. Plot of Hee(8,8) Calculated Using the Synthesis Method. 
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Figure 6.5. Plot of H(8,8) Calculated Using Traditional FE Procedures. 
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Figures 6.4 and 6.5 are the plots of the FRF at element (8.8) calculated using the 
synthesized [F/„j and FE [H\ matrices, respectively The FRF element (8.8) corresponds 
to the lateral motion coordinate 8 of Figure 6.1. Notice both plots are identical and show 
the first five damped natural frequencies. The plots show the magnitude of the response at 
coordinate 8 due to a unit excitation at varying frequency at coordinate 8. As the Irequency 
of excitation approaches the damped natural frequency, the response approaches infinity. 

A comparison of the compute time required for the synthesis versus tradidonal FE 
calculation. The actual computing time and the number of floating point operations (flops) 
for each method is provided: 

FEM direct assembly: time - 25876 sec or 431.3 mins 
FLOPS - 1.49 X 109 

FRF synthesis: time - 1167 sec or 19.45 mins 
FLOPS-517.2 x 106 

This clearly demonstrates that synthesis by FRF is more efficient and well suited for design 
analysis. 

G. EXAMPLE (7): STRESS CALCULATION BY DYNAMIC INDIRECT 

COUPLING 

Consider the structures shown in the following figures. The structure shown in Figure 
7.1 will be directly assembled by the finite element method and the peak bending stress 
frequency response will be calculated in beam element #4 whose location is shown by the 
dashed line A-A. The same structure will be synthesized using the frequency domain 
method and the same stress frequency response will be calculated. The FRF results 
calculated by the synthesis methodology, shown in Figure 7.5, will be compared with that 
calculated by traditional FEM, Figure 7.6. Again, the structure (specifically, its FRF) as 
shown in Figure 7.1 will be obtained by synthesizing structure 1 and the modification. 
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structure 2, as shown in Figure 7.2. Note that a stress frequency response allows the direct 
calculation of stress due to the application of a force or moment. The equation for 
determining synthesized stress is shown as equation (105). 

ja(a))‘ = [//jn)f|/(n)i. dosi 

where the synthesized stress FRF matrix [//^^(O)] reflects the total structure, including 
any modifications or couplings. 



Figure 7.1. Structure Analyzed for Peak Bending Stress 



Figure 7.2. Components of Synthesized Structure 
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Referring to Figure 7.2, structure 1 will be modified by adding the beam, structure 2 at the 
nodal coordinates 4. 5, 6, 10, 11, and 12. The following beam element data will be used: 
Young's Modulus E = .30.0 » 10* psi 
Area moment of inertia I = 0.02083 in^ 

Cross-secaonal area A = 1 in 
Weight density WTD = 0.2832 Ibf/in ' 

Distance from beam center to outer most fiber c = 0.05 in 
For this example, damping was not used, but the methodology is able to handle all 
forms of linear damping as described earlier. The general equation for synthesizing stress 
information by dynamic indirect coupling/modification is 


^ ( 66 ) 

This equation is the first row extracted from the relationship shown as equation (12). Note 
that the synthesis of stresses can be done at the same time as the synthesis of 
displacements. We are here demonstrating just the synthesis of stress information. 

The first step is to generate the [/C] and [M] matrices for structure 1 and the beam 
element shown in Figure 7.2. Next the impedance matrices are generated for each structure 
as [Z(D)] = [^r] Since we are modifying structure 1, the FRF matrix {H\ is only 

calculated for structure 1 of Figure 7.2. The complete process as described here is 
performed over the frequency range of interest. There is basically two sections to this 
process: (1) the partitioning of the [//] matrix into its required sub matrices for the general 
synthesis process and (2) the extraction of the information from the [ H\ matrix and the 
processing of that information to calculate the stress frequency response. 
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Referring to Example (4). we partition [ H] for structure 1 in the same manner. For the 
synthesis of stress frequency response only, only the partitions of [//,J and [w ] are 
required. These partition are shown below. 

‘i t'l 

[«J = I [W„J = c,[//.(r.c,)] 

The internal coordinates "t, "are 1, 2, 3, 7. 8, and 9. The connection coordinates "f," are 
4, 5, 6, 10, 11, and 12. The second part of the process requires all or part of the FRF 
matrix for structure 1, depending on where the external loads are applied. The connection 
coordinates are required for the synthesis process, as always, and internal coordinates are 
required if the stress frequency response which is of interest is associated with an element 
whose nodal coordinates are internal coordinates, i.e. they are not directly associated with 
the synthesis. In this example, all coordinates are used for the stress information. We apply 
a unit load at each coordinate using the following equation, 

|xr=[«]U|. (106) 

where i is an element of the required coordinates. This equation is interpreted as the 
displacements at the structimal system coordinates due to the unit load at the desired 
coordinates of interest, which is the combination of connection coordinates and any internal 
coordinates desired, (xj' is the i'th column of [W] when using unit forces. Using the i'th 
column of [//], we extract the elements corresponding to the beam element that stress 
information is desired for, getting a partitioned form of the i'th column. The complete 
reduced form of the [H] matrix is \H{bc,dc)\,^^^^ , where be are the coordinates of the 
beam of interest and dc is the set of required coordinates we wish to keep. In our example. 
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the set of beam coordinates, "be" is 1. 2, 3. 10, 11. and 12 and the set of required 
coordinates, "dc" is 1,2, 3, 4. 5, 6. 7, 8, 9. 10, 11. and 12. f is a matrix of size 

(6 X 12). Each column of \ is in the global coordinate system and needs to be 

transformed to local coordinates by using the following relation 

(107) 

where the transformation matrix [7] is 

sina 0 0 0 0" 

cosa 0 0 0 0 

0 10 0 0 

0 0 cosa sina 0 

0 0 -sina cosa 0 

0 0 0 0 1 

Now with the H transformed to the local coordinate system, we can get the nodal forces by 
the following relation 

(/'I (108) 

where is the elemental stiffness matrix in the local coordinates system. We are ready to 
solve for the FRF stress, first combining all the column vectors {/') into a nodal force 
matrix [F] and then multiplying it by the moment equation to solve for the FRF peak 
internal bending stress of the beam element 
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The FRF stress equation is 

. ( 109 ) 

The {//„} is a row vector size (1 x 12), since we used all the coordinates as our set ot 
desired coordinates. Noting that jw,,} is determined from equilibrium for the element in 
question, and here we provide the internal bending moment. The derivation is shown using 
Figures 7.3 and 7.4. 



Beam Element 

Figure 7.3. Beam Element for Stress Calculation 


Sectioned Beam Element 

Figure 7.4. Beam Section Cut at the Midpoint 

Consider the beam element in Figure 7.3, The moment is to be calculated at the midpoint of 
the beam. First the beam is cut, Figure 7.4 and the moment at A is solved for: 

M^=M,+V,^ (110) 

The moment equation in vector form is shown as equation (111). 
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Noting that the nodal force of the beam element is 


v;- 



K 


where “A” indicates an axial force, “V” a shear force, and “M." a moment. The internal 
bending stress frequency response component is determined by equation (113). 


w :=-|0 //2 1 0 0 01 


(113) 


//„ is evaluated over all the chosen required coordinates to form {//„}, which in our 
example is a row vector size (1 x 12). If more then one beam element is used for stress 
calculations then [ //„ ] could be of the size; (number of beam elements) x (number of 
desired coordinates). With [//, ] generated, we can now partition it into the required sub¬ 
partitions for synthesis. The partitions required are and [//^J and are 
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= nl\HJnb,cVf\ [H^^\ = nt{Hjnh.i\) H^inb.ci)]' 

The beam element set is indicated by "nb". In this example "nb" could have been 1, 2. 3, 
and 4 since there are four beam elements in the substructure to be modified. Beam element 
4 was chosen as the beam to calculate FRF stress information so "nb" is 4 and the sizes of 
the matrices are (1 x 6) and (1 x 12) respectively. To get the stress mformation for beam 
five, the synthesis method of direct coupling must be used. With the appropriate 
partitioning completed the synthesis can now be performed. [//,J is the modified FRF 
stress relation which is the combination of structure 1 and the added element structure 2 of 
Figure 7.2. The synthesis is performed over the frequency range of interest and plotted in 
figure 7.5. The frequency range for this example is 0.3 to 102 Hz. Figure 7.6 is the 
solution from the traditional FE calculation for direct comparison of the two solutions. Both 
plots are identical. 



Figure 7.5. Plot of Synthesized FRF Stress Element W'„,(l,9) 
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Figure 7.6. Plot of H„{\ .6) Calculated Using Traditional FE Procedures. 

Figures 7.5 and 7.6 are the plots of the FRF stress corresponding to beam element four of 
Figure 7.1. These plots represent the stress amplitude in beam element four due to a unit 
force applied at coordinate 6. Both plots are identical. Note that the frequency of peak 
response is slightly lower than the undamped natural frequency. 

II. EXAMPLE (8): DYNAMIC DIRECT COUPLING USING MODAL 
REPRESENTATION OF ¥Rt 

Consider the structures shown in the following figures. The structure shown in Figure 
8.1 will be directly assembled by the finite element method in order to compare the 
frequency response calculated by traditional FE methods and the solution obtained by 
synthesizing structure 1 and structure 2 as shown in Figure 8.2. This example will show 
three results, the first being the solution from synthesis using the modal representation of 
the frequency response, the second being direct assembly using FE and the modal 
representation and thirdly, direct assembly using FE where frequency response is 
calculated by the inverse of the impedance matrix. 
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Figure 8.1. Stnicaire Analyzed Using Traditional FE 
123 456 123 456 



Figure 8.2. Structures to be Synthesized Using Modal Representation 


Referring to Figure 8.2, structures 1 and 2 will be synthesized by direct coupling using 
connection coordinates 4, 5. and 6 of structure 1 and connection coordinates 1, 2. and 3 of 
structure 2. One internal coordinate will be kept in this synthesis process to show that the 
frequency response for a specific coordinate can be synthesized using just the connection 
coordinates and any internal coordinates that might be of interest This example will use the 
internal coordinate "2" of structure 1 as the coordinate of interest. The information desired 
in this example is the frequency response at coordinate 2 due to a unit harmonic load at 
coordinate 6. Note that when structure 1 and structure 2 are synthesized, the coordinate 
numbering becomes the same as depicted in Figure 8.1. The following beam element data 
was used: 

Young's Modulus E = 30.0 * 10'* psi 
Area moment of inertia I = 0.02083 in" 

Cross-sectional area A = 1 in 
Weight density WTD = 0.2832 Ibf/in* 

Structural proportional damping was not used, but the methodology will handle all 
forms of damping discussed earlier. The frequency response matrix [ //] can be generated 
by two methods. The first is by the relationship 
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[w(0)l = [z(a)f' 


(14) 


where [Z(D)] = [A:]-0'[W] + jC{C] . 


The second method is by matrix modal representation. The relationship is 


[Ma)] = [4>) 


kr 


(42) 


where [4>] is the set of eigenvectors or mode shapes, and the middle term is the diagonal 
matrix of the natural frequencies or eigenvalues less the frequency of interest. This 
relationship can also be expressed in terms of individual elements of the frequency 
response function by equation (114) 



(114) 


which allows the calculation of specific frequency response of interest without having to 
generate the complete FRF. 

The general synthesis equation for dynamic direct coupling is 

= Kl - • (84) 
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This equation is written in terms of all the coordinates. This example concerns the synthesis 
of a single FRF matrix element involving one coordinate. In general, the synthesis process 
requires FRF information for all connection coordinates, and FRF information for any 
internal coordinates of interest. Rewriting the general equation for this specific example, the 
equation becomes 

[^2,]’ (115) 

where the subscript “2” signifies the internal coordinate "2" of structure 1 and the subscript 
“c" signifies the set of connection coordinates of both structures. In this equation, [M] is 
the boolean mapping matrix which is used to establish the connectivity between the two 
substructures for synthesis. The mapping matrix is determined by the connectivity i.e. what 
is connected to what and by imposing the equilibrium and compatibility relations associated 
with each pair of coordinates. We can define the mapping matrix by j/ j = [M]|/,j. 

Where {/.) is a vector of all the connection coordinates of both structures and j/^, | is the 

arbitrarily selected independent subset of the connection coordinates relating to one of the 
substructures. We have selected the connection coordinates of structure 1 as the arbitrary 
subset of connection coordinates. The mapping matrix [M] is a matrix of size (6 x 3) and is 
depicted as: 


[A/] = 
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The FRF matrix [ H] for both substructures 1 and 2 is calculated using the matrix modal 
representation relation discussed earlier. The important point here is that we are not using 
all the coordinates. The coordinates used from substructure 1 are 2, 4, 5, and 6. the first 
coordinate is the coordinate of interest and the rest are connection coordinates which must 
be used. The coordinates used from substructure 2 are just the required connection 
coordinates 1,2. and 3. All six mode shapes for each coordinate are kept for the calculation 
of the FRF matrix. The FRF matrix [//J is calculated by using the appropriate 
partitioning of the modal matrix The diagram of the relation on the next page is 
showing the coordinates kept and the number of modes. 

6 6 2 4 5 6 

2r iT' Ir 1 



The size of [ //, ] is now a (4 x 4) matrix which contains all the necessary information. This 
also shows a significant computational advantage because the size has been reduced from a 
(6 X 6) to a (4 X 4) matrix which requires less computational time to manipulate the matrix. 
[//,] is generated in the same manner. We vrill use all six mode shapes for each coordinate 
kept of substructure 2. The required coordinates are the connection coordinates 1, 2, and 
3. The diagram of the relation is given on the next page showing the coordinates kept and 
the number of modes. 
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Now with h I and hi generated, the two substructures can be synthesized. Referring to the 
example-specific synthesis equation above, the matrices [//,_] and [W j are formed by 
combining hi and h2 by appropriate partitioning. The partitioning is shown on the next 
page. 




_c, A/,(c„q) [0] 

cA [0] H,{c„cA 


Referring to Figure 8.2, "il" is the set of internal coordinates for substructure 1. Since 
coordinate 2 is the only coordinate of interest, the set of internal coordinates is just 
coordinate 2. The set of connection coordinates "cl" consists of 3, 4, and 5 and "c2" 
consists of 1, 2, and 3. With the appropriate partitioning complete, the two structures are 
synthesized together, to form the structure in Figure 8.1, using the case specific form of the 
direct coupling relation 

["2.r (115) 

[//2J is the synthesized FRF by modal representation relation which is the combination of 
both structures. The synthesis is done over the frequency range of interest and plotted in 
Figure 8.3. The frequency range for this example was 0.1 to 80.0 Hz. Figure 8.4 is the 
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solution from traditional FE calculations using the inverse of the impedance matrix to 
calculate the FRF and Figure 8.5 is the solution from traditional FE calculations using the 
modal representation to calculate the FRF. Figures 8.4 and 8.5 are included for direct 
comparison. All three plots are identical. 



Figure 8.3. Plot of Synthesized 1,3) 






Frequency Hz 

Figure 8.4. Plot of H(2,6) from Traditional FE Calculations 
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Figure 8.5. Plot of H(l,4) from Traditional FE Calculations Using Modal Representation 

Figures 8.3, 8.4 and 8.5 are the plots of the FRF at element (1,3), (2,6), and (1,4) 
respectively. These elements corresponds to the lateral motion coordinate 2 of Figure 8.1. 
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A special note here is that the element (2.6) of the FRF generated by FEM is the response at 
coordinate 2. which corresponds to the element (1,3) of the synthesized FRF generated by 
the direct coupling relation using modal representation, and element (1,4) of the FRF 
generated by traditional FE using modal representation. The reason for this is because of 
the partitioning and the coordinates used in the calculation. Care is required here to ensure 
the coordinate of interest is actually being used. The plots show the magnitude of the 
response at coordinate 2 due to a unit excitation at varying frequency at coordinate 6. As the 
frequency of excitation approaches the natural frequency of response, the response 
approaches infinity. 

Figure 8.6 is the plot of the determinant of which shows the natural frequencies of 

the synthesized structure. The frequencies where the plot crosses the axis or equivalently, 
the frequencies for which the det[w„] = 0 correspond to the natural frequencies of the 

synthesized structure. This information is important because it gives the designer a starting 
point on deciding how many modes to keep in the modeling of the system and the 
frequency bandwidth over which to perform the synthesis. Reducing the number of 
retained modes will decrease the computational cost and the computer time required to 
analyze a given design. The number of modes required to accurately model a given 








V. CONCLUSIONS AND RECOMMENDATIONS 


The most important conclusion from this study is that the analysis and re-analysis of 
structural systems is performed most efficiently by working in the frequency domain. It 
was shown in example (6) that synthesis of complex structures was approximately 22 
times faster than the traditional FE methods. The large increase in efficiency means that 
rapid analysis and re analysis of structures can be performed Large scale structural 
analysis can now be looked at in man hours where analysis by traditional FE methods is in 
man days. 

Structural synthesis in the frequency domain provides for an arbitrary order model 
reduction that requires only the coordinates involved in the synthesis and any other 
coordinates that might be of interest. The solution to the reduced model is exact. This is a 
significant point because a 10,000 degree of freedom model can be reduced to a system of 
tens or hundreds of degrees of freedom, significantly improving the computational 
efficiency. 

The frequency response theory allows for the direct synthesis of response information 
of any kind. Using a generalized definition of frequency response, displacement, velocity, 
acceleration, stress, and strain information may be directly synthesized, based on this 
generalization, the theory is an ideal means for doing static and dynamic design re-analysis. 
Static problems are treated as the zero frequency case. 

The frequency domain structural synthesis theory allows for any combination of 
substructure coupling and structural modification to be performed, either simultaneously or 
sequentially. 
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Recommendations to further utilize and demonstrate the theory of frequency domain 
structural synthesis is first, to write computer code that will interface with existing finite 
element codes, for example NASTRANS (MSC, Corp.) or IDEAS (SDRC, Corp.) to 
synthesize substructures in three dimensions using plate, shell, and beam elements with six 
degrees of freedom per node which allows for out of plane analysis. Using the combination 
of plate, shell, and beam elements will more closely approximate actual structures. Second, 
build a scaled prototype of a submarine and equipment cradle and compare the theoretical 
results with experimental results. 
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APPENDIX A MATLAB CODE FOR EXAMPLE ONE 


clear 

% unifineel 

% 

% This program will calculate the eigenvalues,eigenvectors 
% natural frequencies, and frequency response function matrix 
% for a three degree of freedom at each node element. 

% The system is modeled with beam elements that are 
% aligned in the same plane but at any angle (2-D). 


% •. 

% This program works for a 
% coordinates and thus six 


element modeled witl 


% the user must enter the following data to meet the beam conf 
% (E) youngs modulus psi 
% (I) area moment of inertia in*4 
% (WTD) weight density lbf/in"3 
% (A) cross sectional area in"2 
% conductivity [ the node connection mapping ] 

% node coordinates [ cartesian coordinates for each node ) 

% (bb) proportional damping constant 
clear; 

% call the data file 
beam2 
% 

% start the program clock and flops to determine program running 
% time and floating point calculations 

flops(0); 

% calculate the number of beam elements 
a=size(cor,) ; 



% calculate the number of beam elements proportionally damped 
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% 

% calculate the number of nodes. 
ncdes=bii); 

% convert the coordinates in to the correct units (in.) 

% 

% calculate the beam element lengths and beam angles 
% in radians 

IC=con‘i,l); 

;D-ccn(i,2) 

1 i , i) =sqrt ( (coord (ID, 1: -coord i IC, l) i *2+ (coord (ID, 2) - coord (IC, 2 ■ ■ ■ 2 i ,- 
DX(i.=coord(ID,1)-coord IIC,1); 

DY ' i) =coord (ID, 2; -coord (IC, 2) ,- 
if DX(i)>=0 & DY(i)>=0; 

t(1,i)=acos(DX(i)/I (1, i)) ; 
elseif DX(i)<0 & DY(i)>-0; 

t (l,i)=acos(DY(i) /Id, i) ; ♦pi/2,■ 
elseif DX(i)<0 &. DY(i)<=0; 

t (1, i)=acos (abs(DX(i) ) /I (1, id -pi; 
else 

t (1, i)=acos(abs(DY ’ i; ) /I d, i) ) ♦ (3*pi/2) ; 
end; 
end; 

% 

% call trig function 
[c,si=ftrig(t,numel); 

% calculate radius of gyration 
for i=l:numel 

end; 

% 

% create the global matrix which is all zeroes. 

% 

kg=[zeros(nodes*!,nodes*3) 1 ; 
mg=[zeros(nodes*3,nodes*3) 1 ; 

% 

% assemble the elemental matrices to the global matrix. 

% 

for i=l:numel 

[kel,mel] =f elements (l(i) ,WTD(i) ,I(i),E(i),A(i),r(i),c(i),s(i)),- 

w=con(i,2); 

% 

kg(3*v-2:3*v,3*v-2:3*v) = kg (3*v-2 : 3*v, 3*v-2 : 3*v) + keld-.3,i:3); 
kg(3*v-2:3*v,3*w-2:3*w) = kg (3*v-2 : 3*v, 3*w-2 : 3*w) + keld:3,4:6); 
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% Che user muse adjust the global matrix to meet the boundary conditions 


% CO delete rows 
kg( [BCl , :) = [ 1 ; 
kgd( [BC] , :) = [ 1 ; 
mg! [BC] , : ) = [ ) ; 

% CO delete columns 
kg! : , [BCl ) = [ 1 ; 
kgd(: , (BCl ) = [ 1 ; 
mg ( : , (BCl ) = ( 1 ; 

% call the function and calculate eigenvectors and the eigenvalues 
% which are the mode shapes and undamped natural frequency in 
(rad/sec)*2 


(lambda,phil=fgmodes(kg,mg); 

% this now converts the eigenvalues to nat frequency in (rad/sec) and 
% and hertz(1/sec) 

omega = sqrt (lambda) ,- 
freq = omega/( 2 *pi); 
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% ccins-rucc tre £re<iuenry respsnse pip;; over che frequencies 
% of inceresc 

for Con’.ega=2:. 5:37 5 
counc=counc* i 
Z=Kgd■C omega■2 *rag; 


% mis decermines the coordinate of inceresc to plot 

HH(count)=H'5,5); 

end; 

% end the program cloc.'C and flops 

Comega=2:.5:375; 

Frea=Comega/(2*pi!; 

plot (Freq, 20»log(HK) ,1 ,grid 

xlabel(' Frequency Hz ') 

ylabeil'FRF at coordinate of inceresc dB') 

% 

% EMD 


% beam2 

% This is Che data for the full structure II 
% 

% The data will be in the form of 

% - (E) youngs modulus psi 

% - (I) area moment of inertia in"4 

% - (WTD) weight density lbf/in'3 

% - (A) cross sectional area in*2 

% - conductivity [ the node connection mapping 1 

% - node coordinates [ cartesian coords for each node in ft. ] 

% the main program will convert to in.. 

% - (bb) structural proportional damping constant 
% 

E=(30 30 30 ]*ie6; 

I=[.i666 .1666 .1666 l*le-3; 

A=i.2 .2 .2 1 ; 

WTD=;.2832 .2832 .2832 ]; 
bb=.02; 

% nodal connectivity 
con=[i,2; 

2,3; 

3,4] ; 

% nodal damping connectivity 
dcon=[1,2; 

3,4] ; 

% nodal cartesian coordinates 
coord= [0,0; 
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% - node coordinates ■ cartesian coords for each node in ft. 1 
% the main program will convert to in.. 

% - (bb; structural proportional damping constant 
E=(30 l*le6,- 
I=[.1666 ]*ie-3; 

A= [ .2 1 ; 

WTD=[.2832 ]; 
bb=. 02 ; 

% nodal connectivity 
con=[l,2] ; 

% nodal damping connectivity 
dcon=[l,21; 

% nodal cartesian coordinates 
coord=[0,0; 

2,01 ; 

% Boundary conditions 
% structure 1 
%BC=[1 2 3] ; 

% Structure 2 
BC= [4 5 6] ; 

% impedence z 
%BC=[ 1 ; 


% This is example i which demonstrates dynamic indirect coupling. 

% Two structures will be synthesized together by way of the new load 






% Che K and M matrix 5cr earn scruccnrs is saved 
lead ci.mat % kl,ml is scored here 
load cl.mac % k2,m2 is scored here 
irad c3.ir.ac % kz,mz is scored hers 
% 

* we need cr rreace a single FRF macrix representing 
i Doch sufcscrurcures in the farm: 

% 

% :h(c,i! h.r.c): 

% 30 we create arrays containing the DOF numbers of our original 
% models which correspond co cne 'c' and "i" coordinates for 
% each subscruecure. 

% 

% call Che synthesis daca file in now which contains the 
% internal coordinates and connection coordinates for each sub 
% structure. 

% il= internal coords of sub structure l 
% i2= internal coords of sub structure 2 
% cl= connection coords of sub structure 1 
% c2= connection coords of sub structure 2 
% 

FRF_INDIR_DATA 

% 

co=clock,- 
flops(0); 
count=0; 

for Comega=2:.5:375 
% 

% Form Frequency Response Models for Each Substructure 


zi=kl -Comega‘2*ml 
z2=k2-Comega■2*m2; 
z=k3-Comega-2*m3; 



b=size(cl); 
C=size(i2); 
d=size(c2); 



cc=c(2); 
dd=d(2); 

% 

% Remember, we are crying co calculate the following: 
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% 1.- ir.rerna: coords suo scruccure l 
% i2- inrernal coords sufc srruccure 2 
% connection coords sub scructure i 
% c2• connection coords sub structure 2 


't 

% enter the number of unrestrained nodes of the synthesized 


99 






APPENDIX B MATI.AB CODE EOR EXAMPLE TWO 


% unifineell 

% 

% This program will calculate the eigenvalues, eigenvectors 
% natural frequencies, and frequency response function matrix 
% for a three degree of freedom at each node element. 

% The system is modeled with beam elements that are 
% aligned in the same plane but at any angle {2-D) 


% 


% •.• 

% This program works for a beam element modeled with six general 
% coordinates and thus six DOF. 

% ( I -• .•- I ) 

% 

% The user must enter the following data to meet the beam configuration. 
% (E) youngs modulus psi 
% (I) area moment of inertia in*4 
% (WTD) weight density lbf/in*3 
% (A) cross sectional area in*2 
% conductivity [ the node connection mapping ) 

% node coordinates [ cartesian coordinates for each node ] 

% (bb) proportional damping constant 
% call the data file 
minihull_data2A 

% Start the program clock and flops to determine program running 
% time and floating point calculations. 
tO=clock; 
flops(0); 

% 

% calculate the number of beam elements 



% calculate the number of beam elements porportionally damped 
aa=size(dcon); 
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~q ' 3*v-2 : 3*v, 3*v-2 : 3*vi = mg ■ 3*v-2 : 3*v, 3*v -2 : 3 *vi mel(i:3,i:3i; 

rng !3*v-2 : 3*v, 3*w-2 : 3*w; = mg i 3*v-2 : 3*v, 3*w-2 : 3*w; * mel(i:3,4:6); 

mg ( 3*w - 2 : 3*w, 3*v-2 : 3*v) = mg ( 3*w-2 : 3*w, 3*v-2 : 3*v; + mei;4:6,i:3l,- 

mg(3*w- 2 :3*w,3*w-2:3*W) = mg(3*w-2:3*w,3*w-2:3*wi * mel(4:6,4:6); 

% 

% 

% apply structural prop, damping to the k matrix and sec global 
% k matrix to equal damped matrix 
% 

% calculate Che beam element lengths and beam angles in radians 
% for Che damped beams 
for i=irnumeldamp 
IC=dcon(i,1); 

ID^dcon(i,2); 

1(1, i) =sqrc ( (coord (ID, 1) - coord (IC, l) ) *2 + (coord (ID, 2) -coorddC, 2) ) - 2) ; 
DX{i)=coord(ID, i) -coord, ZC, 1) ; 

DY(i)=coord(ID,2)-coord,10,2); 
if DX(i)>=0 & DY(i)>=0; 

C(1,i)=acos(DX(i)/l(i,i)); 
elseif DX(i)<0 s, DY(i)>=0; 

c(1,i)=acos(Dy(i)/I(1,i))+pi/2; 
elseif DX(i)<0 & DY(i)<=0; 

C(1,i)=acos(abs(DX(i))/I(l,i))*pi; 
else 

c(1,i)=acos(abs(DY(i))/l(l,i))+(3*pi/2); 
end; 
end ,- 
% 

% call trig function 
[c,sl=fcrig(c,numel_damp); 

% 

% calculate radius of gyration 
for i=i:numel_damp 
r(i, i)=sqrt(I(i)/A(i)); 


for u=l :numel_daiTip 

[kel] =f element 6 (1 (u) , WTD (u) ,I(u) ,E(u) ,A{u) ,r(u) ,c(u) ,s(u) ) ; 

% 

v=dcon(u,1); 
w=-dcon.(u., 2) ; 

% 

kgd(3*v-2:3*v,3*v-2:3*v) = kgd(3*v-2:3*v,3*v-2:3*v) + j*bb*kel(1:3,1 

kgd(3*v-2:3*v,3*w-2:3*w) = kgd(3*v-2:3*v,3»w-2:3*w) + j*bb»kel(1:3,4 
kgd(3*w-2:3*w,3*v-2:3*v) = kgd(3*w-2:3*w,3»v-2:3*v) + j*bb*kel(4:6,1 
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- ;»bb»> 


.<gd -*>^-2:3*1^, 3*w}*w = 


% apply v.he boundary conditions 

% the osar must adjust cne global x.atrix to meet the boundary conditions 



kgd! ;bc: , : = !; ; 

mg i [BCi , : = .1 ; 

% to delete columns 

kg :, ;Bc: = ; 

kgd> :, IBC) ) = 
mg I : , iBC; > = ; 

% 

% call the function and calculate eigenvectors and the eigenvalues 
i which are the mode shapes and undamped natural frequency in 
%;rad/sec)'2 

Llambda,phi;“fgmodesfkg,mg); 

% convert the eigenvalues to nat frequency in (rad/sec) and 
% and hertz(i/sec) 

% 

omega = sqrt .ambda); 
freq = omega/(2*pi); 

% construct the frequency response plot over the frequencies 
% of interest 

for Comega=.l:.1:22 
count=count+l; 

Z=kgd-Comega"2*mg; 

H=inv(Z); 

% determines the coordinate of intrest to plot 
Kh • count) =;• 8,8) ; 

% end the program clock and flops 
etime(clock,tO),flops 
Comega=.1:.1:22; 

Freq=Comega/(2*pi); 

plot(Freq,20*log(HH)),grid 

xlabel(' Frequency Hz ') 

ylabeK'FRF at coordinate of interest dB ') 



minihull_data2A 
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% degrees of freedom ac a node. 


% 

% data will be in the form of 

% - (E. youngs modulus psi 

% - (I) area moment of inertia in"4 

% - (WTD) weight density lbf/in*3 

% ■ {Ai cross sectional area in"2 

% - conductivity [ the node connection mapping ] 

% - node coordinates [ cartesian coords for each node in ft. ] 
% the main program will convert to in.. 

% - (bb) structural proportional damping constant 


E=[30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 301*le6; 

1=;.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .020831 ;% 1/i2bh*3 b=2, 
h=. 5 

A=[lllllllllllllllllll; 

WTD=[.2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 
.2832 .2832 .2832 .2832 .2832 .2832 .28321; 
bb=0.01; 
con=[l,2; 

2,3; 


8,9; 

9,10; 


determines 


elements have damping 
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% - (A) cross sectional area in*2 
% - conductivity [ the node connection mapping ] 

% - node coordinates [ cartesian coords for each node in ft. 1 
% the main program will convert to in.. 

% - (bb) structural proportional damping constant 
E=[30 30 30 30 30 30 30 30 30 30 30 301*le6; 

I=i .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083];%i/12bh'3 b=2, h=.5 

WTD=[.2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 
.2832); 
bb=0.01; 
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% So we create arrays containing Che DOF numbers of our original 
% models which correspond the the "c" and "i" coordinates for 
% each substructure. 

% 

% Call Che synthesis data file which contains the 

% internal coordinates and connection coordinates for each substructure. 
% 

% il= internal coords of sub structure l 
% i2= internal coords of sub structure 2 
% cl= connection coords of sub structure 1 
% c2= connection coords of sub structure 2 
% 

FRF_Synch_daca2A 

% 

tO=clock;; 
flops(0); 
counc=0; 

for Comega=.1:.3:60 


% Form Frequency Response Models for Each Substructure 


zl-=kl -Comega‘2*ml ; 
z2=k2-Comega*2*m2; 



aa=a(2); 



108 







se tnacri-es connain zhs FRF data for both subscrucnu 
or CO coupling, i.e cne pre-synchesis FRF data. 


Coordinate Partitioning 


id up uncoupled FRF matrix and sub-partitions: 

,ni'li,ii) zeros:aa,cc) hlui,cij zeros (aa,dd: ; 
zerosicc.aai h2 i2,i2) zeros(cc,bb) h2ii2,c2'; 
n-ioi.ii. zeros bb,CO hi(ci,cii zeros i bb, dd) ; 
zeros dd,aa) h2’C2,i2 zeros(dd,bb) h2(c2,c2;i; 


Ihl!il,cl; zerosaa.dd.i; 
zerosiCC,bbi h2-12,02 
hi(cl,ci; zeros>bb,ddi; 
zeros{dd,bb) h2ic2,c2)|; 


;hi'’cl,cl) zeros (bb,dd) ; 
zeros(dd,bb) h2(c2,c2)); 

[hl(cl,il) zeros(bb,cc) hi(cl,cl) zeros(bb,dd) 
zeros(dd,aa) h2{c2,i2) zeros(dd,bb) h2(c2,c2)l; 


% We can now perform the synthesis: 
hccr=M' ♦ hcc * M; 

heescar = hee - hec * M * inv( hccr ) » M' * hce; 


% remove the redundenc information 
heescar = heescar(i:nodes*3,l:nodes*3); 
% loo)< at the coordinate of interest 
HH(count)=heestar(8,8); 


ecime(clock,cO),flops 
C:omega=. 1: . 3 :60 ; 

Freq=Comega/(2*pi) ; 

plot(Freg,20*log(HH)),grid 

xlabel(’ Frequency Hz ') 

ylabeK'FRF at coordinate of interest dB') 


% FRF_Synch_daca2A 

% This is Che data file for the synthesis program. 
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Al’PFNDIX C MATLAB CODK FOR EXAMPI.E TffREE 



% unifinee. 

% 

% This program will caloulace The eigenvalues,eigenvectors 
% natural frequencies, and frequency response function matrix 
% for a three degree cf freedom at each node element. 

% The system is modeled with beam elements that are 
% aligned in “f'e same plane but at any angle (2-Dt . 


% This program works for a beam element modeled with six general 
% coordinates and thus six DOF. 


% the user must enter the following data to meet the beam configurat 
% (E) youngs modulus psi 
% (I) area moment of inertia in'4 
% (WTD) weight density lbf/in'3 
% (A) cross sectional area in'2 
% conductivity [ the node connection mapping ] 

% node coordinates I cartesian coordinates for each node 1 
% (bb) proportional damping constant 

% call the data file 
minlhull_data3 

% start the program clock and flops to determine program running 

% time and floating point calculations 

tO=clock; 

% calculate the number of beam elements 

numel^a(l); 

% 

% calculate the number of beam elements porportionally damped 
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% calculate the number of ncdes. 
b=size(cocrd); 


% 

% convert the coordinates in to the correct units (in.) 

% calculate the beam element lengths and beam angles in radians 
for i=i:numel 



1 (1, =sqrt((coord(ID,1)-coord(IC,1))*2+(coord(ID,2)-coord(IC,2)) -2); 
DX(i)=coord(ID,1)-coord(IC,1); 

DY(i:=coord(ID,2) -coord(IC,2) ; 
if DX(i) >=0 & DY(i) >=0; 

t (1, i)=acos(DX(i) /Id, i) ) ,- 
elseif DX(i)<0 & DY(i)>=0; 

t(1,i)=acos(DY(i)/I(1,i))*pi/2; 
elseif DX(i)<0 & DY(i)<=0; 

t d, i)=acos (abs(DX(i) ) / I (1, i) ) tpi ; 
else 

t (1, i) =acos (abs (DY(i))/l(l,i)) + (3*pi/2) ; 
end; 

% 

% call trig function 
[c,sl=fcrig(t,numel); 

% 

% calculate radius of gyration 
for i=l:numel 
r(1,i)=sqrt(I(i)/A(i)); 
end; 

% 

% create the global matrix which is all zeroes. 

kg=[zeros(nodes»3,nodes»3)]; 
mg=[zeros(nodes*3,nodes*3) 1 ; 

% 

% assembel the elemental matricies to the global matrix, 
for i=l:numel 

[kel, mel) =f element 6 (1 ( i) ,vgTD (i) ,I(i) ,E(i) ,A(i) ,r(i) ,c(i) ,s(i) ) ; 



kg(3*v-2:3*v,3*v-2:3*v) = kg(3*v-2:3»v,3*v-2:3*v) + kel(l:3,l:3) 
kg(3*v-2:3*v,3*w-2:3*w) - kg(3*v-2:3*v,3»w-2:3*w) + kel(l:3,4:6) 
kg(3*w-2:3*w,3*v-2:3*v) = kg(3*w-2:3*w,3*v-2:3*v) + kel(4:6,l:3) 
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3,1:3'; 




3*v = mg , 3‘v - 2 : 3‘v, 3*v-2 : 5*v; * me* 1. 

3*wi = mg,3*v-2:3*v,3*w-2:3*wi * mel' 1 : 3, 4 :6 , ; 
3*V' = mg ; 3*w-2 : 3*w, 3*v-2 : 3*v) - irel 4:6,1:3); 

3*w) = mg{3*w-2:3*w,3*w-2:3*w) raei 14 :6,4 :6 ) ; 


% 

% apply structural prop, damping to the k matrix and se' gioCal 
% k matrix to equal damped matrix 
% 

% 

% talculate the beam element lengths and beam angles in radians 
% for the damped beams 

IC=dcon(1,1.; 

:h-dcon 11,2i; 

1; i, 1 • =sqrt (coord i ID, ii - coord dC, i) i '2* (coord i ID, 2) ■ coord (IC, 2, ) * 2 
DX11;=coordI ID,1)- coord11C,1! ; 

DY i, =coord(ID,2)-coordiIC,2); 
if DX(i)>=0 & DY(il>=0; 

t 1, i.'=acos(DX(i) /Id, i, ) ; 
elseif DXdl<0 & DY(ii>=0; 

t (1, ii =acos (DY( 1) / I ■; i, i) ) -pi/2,- 
elseif DXdXO & DY(i)<=0; 

t(1,i)=acos(abs(DX(i))/I(1,1))-pi; 


t(l,i; --acos(abs(DY(i))/l(l,i)) + (3*pi/2) ; 
end; 


% call trig function 
;c,s]=ftrig(c,numel_damp) ,- 
% 

% calculate radius of gyration 
for i=l:numel_damp 
r(1,1)=sqrt(I(ii/A(i)) ; 
end; 


kgd=kg; 

for u-l:numel_damp 

tell =felements(l(u),WTD(u),I(u),E(u),A(u),r(u),c(u),s(u)); 

v=dcon(u,1); 
w=dcon(u,2); 

% 

kgd(3*v-2:3*v,3*v-2:3*v) = kgd(3*v-2:3*v,3*v-2:3*v) - j*bb*kel(i:3, 

kgd{3*v-2:3*v,3*w-2:3*w) =■ kgd ( 3*v-2 : 3*v, 3 *w-2 : 3*w) - j *bb*kel (1: 3 , 
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% call the function and calculate eigenvectors and the eigenvalues 
% which are the mode shapes and undamped natural frequency in 
(rad/sec)'2 

' lambda,phi] =fginodes (kg,mg) ; 


end; 

% end Che program clock ai 
ecime(clock,CO),flops 
Comega=.l:.3:40; 
Freq=Comega/(2*pi); 
plot (Freq, 20*log (HH) ) ,gri( 
xlabeK' Frequency Hz ') 
ylabel('FRF at coordinate 
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= daia will be in the form of 
(E. youngs modulus psi 
I) area moment of inertia in"' 
WTD, weight density lbf/in*3 


% ■ node coordinates 


;esian coords for each node in ft. ; 
main program will convert to in.. 

% - bb) structural proportional damping constant 
i 

E=i 31 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30J*Le6; 

I=i .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083);%i/l2bh-3 b=2, 


.2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 


)32 .2832 .2832 .2832 .2832); 


deteremines what eleri 


7,8; 

8,9; 

9,10; 

10 , 11 ; 

11 , 12 ; 

12 , 1 ; 

11,13; 
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% This is the data file for the beam modification. 

% The data will be in Che form of 
% - (E) youngs modulus psi 


% - (WTD) weight density lbf/in*3 
% - (A) cross sectional area in*2 
% - conductivity [ the node connection mapping 1 
% - node coordinates [ cartesian coords for each nods 
% the main program will convert 

% - (bb) structural proportional damping constant 
% 

E=[30 ]*le6; 

I=[. 02083 ] ; 


bb=0.01; 

% this determines what elements have damping 
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% 

% rr.;3 is -r.s iaia file for ohe T.ain sur.soroo'ure. 

% The data will be in the form nf 

i ■ 'E' youngs modulus psi 

A - 111 area moment of inertia in*4 

% - iWTD' weight density ibf/in*3 

% - Ai rross sectional area in*2 

% - conductivity [ the node connection mapping ) 

% - node coordinates [ cartesian coords for each node in ft. ] 

% the main program will convert to in.. 

% - bb) structural proportional damping constant 

E-:30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30|*ie6; 

1=1.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .020831;%!/12bh*3 b=2, 

WTD=[.2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 

.2832 .2852 .2832 .2832 .2832 .2832 .2832]; 

bb=0.01; 



is determines 


elements 


damping 
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ihecl 



hcc = [h(cc,cc)); 

% 

hce = [h(cc,ic) h(cc,cc)l; 

% this is for adding a component 
%heescar = hee - hec * inv(inv(2) <■ hcc) * hce; 
% 

% this is for removing a component 

heestar = hee - hec • inv(hcc - inv(z)) * hce; 

HH(count)=heestar(11,11); 
end; 

etime (cloclt, to) , flops 
Comega=.1;.3:40; 

Freq=Comega/(2*pi) ; 

plot(Freq,20*log(HH) ) ,grid 

xlabel(' Frequency Hz ') 

ylabeK'FRF at coordinate of interest dB ') 


% MOD_DAT3 

% This is the data file for the modification program, 
% the following data will be provided by this file 
% 

% ic- internal coords of synthesized structure 
% cc- connection coords of synthesized structure 
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ic=[i 2 3 4 5 6 7 8 9 13 14 15 L6 17 18 .3 20 21 22 23 24 25 26 27 28 23 
30 31 32 33 34 35 36 37 38 39 40 41 42 ]; 

CC=flO 11 12 43 44 451 ; 
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APPENDIX D MATLAB CODE EOR EXAMPLE POUR 


% This program will caicui.ace the eigenvalues, eigenvectors 












:-g 3*v-2 : 3*v, 3*v-2 : 3*v = mg i 3*v-2 : 3‘v, 3*v- 2 : 3*v, * mel,i:3,J.:3); 

mg ' 3*v-2 : 3*v, 3*w-2 : 3*w, = mg i 3*v-_ : 3*v, 3*w2 : 3*w.' • mei i i : 3,4 : 6 ! ; 

mg 3*w-2:3*w,3*v-2:3*v, = mg(3*w-2:3*w,3*v-2:3*v - mel(4:6,1:3j ; 

mg 3*w 2:3*w,3*w-2:3*w) = mg(3*w-2:3*w,3*w-2:3‘wl - mei14:6,4:6); 

end 
% 

% apply scruccural prop, damping zo the k .macrix and sec global 
% k macrix co equal damped macrix 
% 

% caicuiace che beam eiemenc lengchs and beam angles in radians 
% for che damped beams 
for i=l:numel__damp 
12=dcor. ’1,1) ; 

1 1, i; =sqrc ((coord(ID, l; - coord: 22, i) i ■2* (coord(ID, 2i -coorddC, 2) ; *2) ; 
DX ( ij =coord .ID, 1) - coorddC, II ; 

DY ( i; ^coord (ID, 2) -coorddC, 2; ; 
if DX(i,>=0 i DY(i)>=0; 

c (1, i)=acos iDXdj /I (1, i) ) ; 
elseif DX(i)<0 i DY(i)>=0; 

c(1,iI-acos(DY(i)/l(l,i))^pi/2; 
elseif DX(i)<0 & DY(i)<=0,- 

c (1, i)=acos (abs (DX( i) ) / I (1, i) ) +pi,- 
else 

c(1,i)=acos(abs(DY(i))/l(l,i))+(3*pi/2); 
end; 
end; 

% 

% call crig funccion 

[c, s] =f Crig (c ,numel_darap) ,- 


% caicuiace radius of gyracion 
for i=1:nume 1 _dan:ip 
r(l,i)=sqrc(I(i)/A(i)); 

% 


kgd=kg; 

for u=l:numel_damp 
[kell=felemenc6(l(u) ,wtd(u) 
% 

v=dcon(u,1); 
w=dcon(u,2); 

% 

kgd(3*v-2:3*v,3*v-2:3*v) = 
kgd(3*v-2:3*v,3*w-2:3*w) = 
kgd(3*w-2:3*w,3*v-2:3*v) = 
kgd(3*w-2:3*w,3*w-2:3*w) = 


,l(u; ,E(u) ,A(U) ,r(u) ,c(u) 


kgd(3*v-2:3*v,3*v-2:3*v) 
kgd(3*v-2:3*v,3*w-2:3*w) 
kgd(3*w-2:3*w,3*v-2:3*v) 
kgd(3*w-2:3*w,3*w-2:3*w) 


+ j*bb*kel(l:3, 
♦ j*bb*kel(i:3, 
+ j*bb*kel(4:6, 
+ j*bb*kel(4:6, 


3) ; 
6) ; 
3) ; 
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boundary cc 


% apply the boundary conditions 
% 

% to delete rows 
kg( [BC! , ; ) = fi ; 

tig; iBC] , : . = ; 

% to delete columns 
kg( : , [BC] ) = [) ; 
kgd! : , [BC) ) = [) ; 

!Tig( : , [BC] ) = [) ; 


% call the function and calculate eigenvectors and the eigenvalues 
% which are the mode shapes and undamped natural frequency in 
% (rad/sec)"2 

[lambda,phi) =fgmodes (kg,mg) ,- 

% this now converts the eigenvalues to nat frequency in (rad/sec) and 
% hertz(1/sec) 

% 

omega = sqrt(lambda); 
freq = omega/{2*pi); 

% 

% constuct the frequency response plot over the frecpuencies 
% of interest 
count=0; 

for Comega=.i:.3:53 

Z=kgd-Comega■2 *mg; 

H=inv{z); 

% this determines the coordinate of intrest to plot 
HH (count) =H ( 8,8) ,- 

% end Che program clock and flops 
ecime(clock,CO),flops 
Comega=.l:.3:53; 

Freq=Comega/(2*pi) ; 

plot{Freq,20*log(HH)),grid 

xlabel(' Frequency Hz ') 

ylabeK’FRF at coordinate of interest dB') 


% minihull_data4 


% This is the data file for the finite element program. 
% 

% Che data will be in Che form of 
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% (E/ youngs modulus psi 
% (I; area moment of inert;ia in'4 
% Ai cross sectional area in"2 
% conductivity [ the node connection mapping j 
% node coordinates i cartesian coordinates for each node i 
% ;bb; proportional damping constant 
% 

E-;30 3i 30 30 30 31 30 30 30 30 30 3C 30 30 30!*ie6; 

I>(.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .020831;%!/i2bh-3 c=2, 
h-.5 

WTD=,.2832 .2332 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 
.2832 .2832 .2832 .2832 .2832 .2832 .2832); 
bb=0.01; 

COn=[l,2; 

2,3; 


elements have damping 


9,10; 
10 , 11 ; 
11 , 12 ; 
12 , 1 ; 
11 , 13 ; 
13, 14; 
14,15; 

13,12; 
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youngs modul'- 


psl 


S (A) cross sectional area in"2 

i conductivity [ the node connection mapping 1 

S node cocrdinates i cartesian coordinates for each node ; 

i ihb) prtpcrtional damping constant 

ii 

E-l 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30i*le6; 

[=[ .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083!;%l/i2Pn*3 


2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 . 


2832 .2832 .2832 .2832 .2832]; 


; elements have damping 


9,10; 

11 , 12 ; 
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for Comega=.l:.3:53 
counc=count i ; 


z=k2•Comega '2*x2; 



hicc.cc!1; 


i is for adding a component 


HH i count)=heescar(8,8); 


ecime(clock,CO),flops 
Comega=.1:.3:53; 

Freq=Comega/(2*pi); 

plot(Freq,20*log(HH)> .grid 

xlabel(' Frequency Hz ') 

ylabell'FRF at coordinate of interest dB ') 


% MOD_DAT4 

% 

% This is the data file for the synthesis program 
% Che following data will be provided by this file 
% 

% ic- internal coords of synthesized structure 
% cc- connection coords of synthesized structure 

ic=[l 2 3 4 5 6 7 8 9 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 
30 31 32 33 34 35 36 37 38 39 40 41 42 ]; 

CC=(10 11 12 43 44 451; 

% 

% enter the number of unrestrained nodes of Che synthesized 
% structure 
nodes = 15; 


129 





,\PPENi:)IX E MATLAB CODE FOR EXAMPLE FIVE 











porp.i 


:ulace rhe number of I 


nume_larnp=aa (1' ; 

% 

% calculate che number of nodes. 
b=3ize;coord ; 

% 

% convert the coordinates in to the correct uni 
coord-coord*12; 




% calculate the beam element lengths and beam angles in radians 
for i=l:numel 
IC-cor. (i,i); 

ID=conil,2i; 

1'1,i;-sqrt((coord(ID,1)-coord(IC,1))*2*(coord(ID,2)- coord(IC,2))"2‘; 
DX(i)-coordiID, i) -coorddC, 1) ; 

DY(if-coord(ID,2)-coordiIC,2,; 
if DX(i)>-0 a DY(i)>=l; 

t (1,i;=acos(DX(i, /I(1, i) ) ; 
elseif DX(i;<o a DY(i)>=:; 

t(l, i)=acos(DY(i) /I 1, i) ) *-pi/2; 
elseif DX(i)<0 a DY(i)<=0; 

t(1,i)=acos(abs(DX(i))/l(i,i))*pi; 

c(l,i)=acos(ab3(DY(i)l /l(i,i))*(3*pi/2); 
end; 

% 

% call trig function 
[c,sj=ftrig(t,numel) ; 

% 

% calculate radius of gyration 
for i=i:numel 
r(1, i)=sqrt(I(i)/A(i)); 
end; 

% create the global matrix which is all zeroes. 

kg=(zeros(node3*3,nodes*3)1; 
mg=[zeros(nodes*3,nodes*3)]; 

% assemble the elemental matricies to the global matrix. 

% 

for i=i:numel 

[kel, men =f elements (1 (i) ,VJTD(i) ,I(i),E(i),A(i),r(i),c(i),s(i)); 
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;3*v, 


3*v, 


kg ( 3*'v- 2 
kg(3*v-2 
kg(3*w-2 
kg(3*w-2 


3*w,3*v- 


2:3*v: 
2:3*w) 
2:3*vl 
2:3*wj 


kg(3*v- 


2:3*v,3*v- 
2:3*v,3*w- 
2:3*w,3*v- 
2;3*w,3*w- 


; 3»w) 


kal a:3, i : 3 j 
kel(i:3,4:6) 
kel(4:6,1:3 j 
kel(4:6,4:6) 


mg(3*v-2 
mg( 3*v-2: 

mg(3*w-2 


= mg(3 *v 
= mg(3 *v 
= mg(3*w- 
= mg(3 


mel(1:3,1:3) 

mel(4:6,1:3; 
mel(4:6,4:6) 


% Tr.is section will apply structural prep, damping to Che k 
% matrix and sec global k matrix to equal damped matrix 

% calculate the beam element lengths and beam angles in radians 
% for Che damped beams 
for i=l :numel_dainp 
IC=dcon(i,li; 

ID=dcon(i,2 ; 

1(1, i)=sgrc'(coord(ID,l)-coord(IC,1))"2+(coord(ID,2)-coord(IC,2) 
DX(i)=coord(ID,1)-coord(IC,1); 

DY(i)=coord(lD,2)-coord(lC,2); 
if DX(i)>=0 & DY(i)>=0; 

c (1, i)=acos (DX(i) /I (1, i) ) ; 
elseif DX(i)<0 & DY(i)>=0; 

t(1,i)=acos(DY(i)/I(l,i))*pi/2; 
elseif DX(i)<0 & DY(i)<=0; 

t(1,i)=acos(abs(DX(i))/I(l,i))+pi; 
else 

t (1, i)=acos(abs(DY(i) ) /1( 1, i) ) + (3*pi/2) ,• 
end,• 
end; 

% 

% call trig function 
[c,sj=fcrig(c,numel_damp); 

% calculate radius of gyration 
for i=l: numel_danTp 
r(1,i)=sgrc(I(i)/A(i)); 

for u=l :numel_dainp 

[kel] =f element6 (1 (u) ,WTD (u) , I (u) , E (u) , A (u) , r (u) , c (u) , s (u) ) ; 

v=dcon(u,1); 

■w=dcon (u, 2) ; 

kgd(3*v-2:3*v,3*v-2:3*v) = kgd(3*v-2:3*v,3*v-2:3*v) + j*bb*kel(l 
kgd(3*v-2:3*v,3*w-2:3*w) = kgd(3*v-2:3*v,3*w-2:3*w) + j*bb*kel(i; 
kgd(3*w-2:3*w,3*v-2:3*v) = kgd(3*w-2:3*w,3»v-2:3*v) + j*bb*kel(4 
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Z:3*w, 


i-2:3*w,3’ 




end 

% 

% nhis section will connecc a spring-damper sysrem Co che g.ocal 
% sniffness matrix, the spring-damper system is made up of a set of 
three springs and dampers chat correspond to che degrees of freedom 
% at a ntde. It attaches to che global stiffness matrix based cn c.ne 
% spring damper connectivity. 

kgds=kgd; 

numspg = di1i; 
counc=0; 

for Comega=. j.: .2:50 
count=councti; 
for j=l:numspg; 

Ikdsprg!=fsprngdamp'k;j) , Comega,B(j ' ); 

x=sdcon (■ j , 1 i ; 
y=sdcon(j,2i; 

kgdsi3*x 2:3*x,3*x-2:3»x) = kgd(3*x-2:3*x,3*x-2:3*x) * kdsprg(i:3,l:3); 
kgds (3’»x-2: 3»x, 3*y-2 : 3*y) = kgd (3*x-2 : 3*x, 3»y-2 : 3*y) + kdsprg(l;3,4:6.,- 
kgds(3*y-2:3<'y,3*x-2:3*x) = kgd(3*y-2:3*y,3*x-2:3*x) t kdsprg(4:6,i:3); 
kgds(3*y-2:3*y,3*y-2:3*y) = kgd(3*y-2:3*y,3*y-2:3*y! + kdsprg(4:6,4:6); 
end 
% 

% apply che boundary conditions 

% che user must adjust che global matrix to meet che boundary conditions 

% to delete rows 
kgds( [BCl , :) = [] ; 
mg((BC!,:) = [1; 

% to delete columns 
kgds ( : , [BC] ) = [1 ; 
mg(: , (BC) ) = (1 ; 

% 

Z=kgds-Comega•2 *mg; 

H=inv(Z); 

% 

HH(count)=H(8,8) ; 
end; 

% 

[lambda,phi)=fgmodes(kgds,mg); 

% this now converts che eigenvalues to nac frequency in (rad/sec) and 
% hertz(1/sec) 

omega = sqrc(lambda); 
freq = omega/(2»pi); 
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ecime:clock,cO),flops 


% 

rreq=Comega/(2*pi); 

plot(Freq,20*log(HH)) .grid 

xlabel(' Frequency Hz ') 

ylabeli'FRF ac cocrdinace of interest dB 


m i n ihu 11 sp r da[n_da c a 5 


This is the data for the finite element program. 

the data will be in the form of 
(E) youngs modulus psi 
(i; area moment of inertia in*4 

(A) cross sectional area in*2 
conductivity [ the node connection mapping ] 

node coordinates [ cartesian coordinates for each node 1 
(bb) proportional damping constant 
(k) spring constant Ibs/in 

(B) proportional viscous damping c 


E=[30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 301*le6; 

I=[.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 . 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02 0831 ;% 1/12bh*3 


A=[l 


. 1 J 


1 1 J 


- 11 ; 


.2832 . 


WTD=[.2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 , 

.2832 .2832 .2832 .2832 .2832 .2832 .28321; 

B=(.01 .01 .01 .01]; % porp damp 2% 

k=[25 25 25 251; % Ibs/in 

bb=0.02; % Structural porp damp 2% 
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% FRF_:;OS?RE'AM_5 

% 

% Lead data from running ur.ifineel program for each substructure. 
% The K and M matrix for eacn structure is saved, 
load exsa.mat % kl ml are stored here 

load exSb.mat % k2 m2 are stored here 

% We need to create a single FRF matrix representing 


% (hee) = ( h(i,i) h(i,c) 1 

% [ h(c, ii h(c,c) ) 

% 

% So we create arrays containing the DOF numbers of our original 
% models which correspond the the *c" and "i" coordinates for 
% each substructure. 

% 

% call the synthesis data file in now which contains the 
% internal coordinates and connection coordinates for each sub 



% il= internal coords of sub structure 1 
% 12= internal coords of sub structure 2 
% cl= connection coords of sub structure 1 
% c2= connection coords of sub structure 2 


FRF_Indsprdam_data5 

% 

to = clock; 
flops(0); 
count=0; 

for Comega=.1:.2:50 

% Form Frequency Response Models for Bach Substructure 


zl=ki-Comega"2*ml; 
z2=k2-Comega"2*m2; 
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Q=sizeic2i; 
% 


dd=d 2 


% Rerrember, we are crying "o calculate the fol*cwing: 

% ;-.ee* = hee - hec » M * inv(zr * hccr ) * M' » hce 

% hccr = M’ * hcc * M 

% zr = pinv(M) ♦ z • pinviM') which is just idencicy matrix size 
% 3 times the number of spring-damp systems 

% So we need to assemble [hee], (hec], [hce] and [hcc] using the 
% the coordinate secs we just defined. 

% These matrices contain the FRF data for both substructures 
% prior to coupling, i.e the pre-synthesis FRF data. 

% Coordinate Partitioning 


Build up uncoupled FRF matrix and sub-particions: 


[hi(il,il) zeros(aa,cc) 
zeros(cc,aa) h2(i2,i2) 
hl(cl,il) zeros(bb,cc) 
zeros(dd,aa) h2(c2,i2) 


hi(il,cl) 
zeros(cc,bb) 


zeros(dd,bb) 


zeros(aa,dd); 
h2(i2,c2); 
zeros(bb,dd); 
h2(C2,C2)]; 


hec - [hl(il,cl) zeros (aa,dd) 
zeros(cc,bb) h2(i2,c2); 
hi(cl,cl) zeros(bb,dd); 
zeros(dd,bb) h2{c2,c2)]; 

% 

hcc = [hl(cl,cl) zeros(bb,dd); 

zeros(dd,bb) h2{c2,c2)]; 

% 

hce= [hi(cl,il) zeros(bb,cc) hl(cl,cl) zeros(bb,dd); 
zeros(dd,aa) h2(c2,i2) zeros(dd,bb) h2(c2,c2)]; 

% 

% We can now perform the synthesis: 
zr = ()c + j*Comega*B»lt) * eye(12); 

heescar = hee - hec * M * inv(inv(zr) + hccr ) ♦ M' ♦ hce; 

HH(count)=heescar(8,8); 

end; 

ecime (cloclt, C0> , flops 
% 

Comega=.l:.2:50; 

Freq=Comega/(2‘pi); 
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APPENDIX F MATLAB CODE FOR EXAMPLE SIX 



% f inesprdanip? 


% This program will calculate Che eigenvalues, eigenvectors 
% natural frequencies, and frequency response function matrix 
% for a three degree of freedom at each node element. 

% The system is modeled with beam elements chat are 
% aligned in the same plane but at any angle (2-D). 


% 


% 


% This program works for a beam element modeled with six general 
% coordinates and thus six DOF. 

% ( I -•-•- I ) 


% the user must enter the following data to meet the beam configu 
% (E) youngs modulus psi 
% (I) area moment of inertia in'4 
% (WTD) weight density lbf/in"3 
% (A) cross sectional area in'2 
% conductivity [ Che node connection mapping 1 
% node coordinates [ cartesian coordinates for each node 1 
% spring-damper conductivity 
% (k) spring constant 

% (q) viscous frequency dependent damping coefficient 
% 

% call the data file 
minihullsprdam_d7 

% calculate the number of beam elements 

C0=clock; 

flops(0); 

numel=a(l); 

% 

% calculate Che number of nodes. 
b=slze(coord); 
nodes=b(l) ; 
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% convert che coordinaces in t: : che correct units ;in.) 
ccord=cocrd*12; 

% 

% calculate the beam element lengths and beam angles in radians 
for i=l:numel 
IC=con(i,1); 

ID=con(i,2;; 

1,1,1 =sdr- coord!ID,1)- coordllC,1))*2+(coord(ID,2;-coord(IC,2)i'2); 
DX(i)=coord(lD, 1) - coord (IC, 1) ; 

DY(i)=coord(ID,2)-coord(IC,2); 
if DX(i)>=0 & DY(i)>=0; 

11i, i)=acos(DX(i)/I(I,i)); 
elseif DX(i)<0 i DY(i)>=0; 

t(1,i)=acos(DY(i)/I(L,i))+pi/2; 
elseif DX(i)<0 & DY(i)<=0; 

t(1, i)=acos(abs(DX(i))/I(l, i)> +pi; 
else 

t (1, i)=acos (abs(DY(i) )/l(l,i))^-(3*pi/2) ; 
end; 
end; 

% 

% call trig function 


% calculate radius of gyration 
for i=L:numel 
r(1,i)=sqrt(I(i)/A(i)); 
end; 


% create the global r 


< which is all zeroes. 


% assemble the elemental matricies to the global matrix. 


I ,A(i) ,r(i) ,c(i) ,s(i) ) ; 


mg(3♦V■ 
mg (3*v 
mg(3*W' 
mg(3*w- 


3*v) = kg(3*v- 
3»w) = kg(3»v 
3*v) = kg(3*w 
3*w) = kg(3*W' 

3*v) = mg (3*V' 
3*w) = mg(3*V' 
3*v) = mg(3*w 
3*w) = mg(3»w 


3*v) ^ keld 
3*w) + keld 
3*v) kel(4 
3*w) + kel(4 

3*v) meld 
3*w) + meld 
3*v) + mel(4 
3*w) + mel(4 
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% This section will connect a spring-damper system the global 
% stiffness matrix. The spring-damper system is made up of a set of 
% three springs and dampers that correspond to the degrees of freedom 
% at node. It attaches to the global stiffness matrix based on the 
% the spring-damper connectivity. 

% 

xgds=kg; 
d = size(scorn; 
numspg - dl1); 
co'unc=C ; 

for j=l:numspg; 

ikdsprg)=fsprngdampCik(j),Comega,g(jy ); 


y=scon:j,2); 

% 

kgds(3*x-2:3*x,3*x-2:3*x) 
kgds(3*x-2:3*x,3*y-2:3*y) 
kgds(3*y-2:3*y,3*x-2:3*x) 
kgds(3*y-2:3*y,3*y-2:3*y) 
end 
% 

% apply the boundary condi 
% the user must adjust the 
% 

% to delete rows 
kgds ( (BCl , :) = !] ; 
mg( [BCl , :) = [] ; 

% to delete columns 
kgds ( : , [BCl ) = [) ; 
mg{: , [BCl ) = [] ; 


kg(3*x-2:3*x, 

kg'3*x-2:3*x, 

kg(3*y-2:3*y, 

kgl3*y-2:3*y, 


3*x-2:3*x) 

3*y-2:3*y) 

3*x-2:3*xj 

3»y-2:3*y) 


kdsprg _:3,1:3) 
kdsprg11:3,4:6 
kdsprg(4:6,i:3 
kdsprg(4:6,4:6! 


to meet the boundary conditions 


% call the function and calculate eigenvectors and the eigenvalues 
% which are the mode shapes and natural frequency in (rad/sec)*2 

lambda,phil=fgmode3(kgds,mg); 

% now convert the eigenvalues to nat frequency in (rad/sec) and 
% hertz(1/sec) 

% 

omega - sqrt(lambda); 
freq = omega/(2‘pi); 


Z=kgds-Comega■2‘mg; 
H=inv(Z); 

% 

HH(count)=H(8,8); 
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% - node coordinates [ cartesian coords for each node in ft. 1 
% the main program will convert to in.. 

% - (q) frequency dependent viscous damping coefficient 
% - (k) spring constant 

E=[30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30]*ie6; 

I=[.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083 
.02083 .02083 .02083 .02083 .02083 .02083 .02083 .02083] ;%i/12bh"3 b=2, 
h=. 5 

A=[iliiillllllllliilll; 

WTD=[.2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 .2832 
.2832 .2832 .2832 .2832 .2832 .2832 .2832); 
q=I.1 .1 .i .i!; % damping coefficient 

k=[25 25 25 25]; % Ibs/in 

% 

COn=[l,2; 

2,3; 


5,6, 
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% 

E=;30 33 30 30 30 30j *196; 

:=[ .:2383 .02083 .02083 .02083 .02083 .02083;;%i/12bh-3 b=2, h=.5 
A=[ 1 1 i 1 i 11 ; 

WTD= .2832 .2852 .2832 .2832 .2832 .2832]; 


COn=i1,2; 



3,0; 
16,0. ; 


3C=[! ; 


clg; 

% 

% FRF_INDSPRDAM_7 

% 

% Load data from running unifineel program for each substructure. 
% The K and M matrix for each substructure is saved. 

% 

load ex7a.mat % kl ml are scored here 
load ex7b.mat % k2 m2 are scored here 
% 

% We need to create a single FRF matrix representing 
% both substructures in the form: 

% [heel = ( h(i,i) h(i,c) 1 

% [ h(c,i) h(c,c) J 


% So we create arrays containing the DOF numbers of our original 
% models which correspond the the "c" and 'i* coordinates for 
% each substructure. 

% call the synthesis data file in now which contains the 
% internal coordinates and connection coordinates for each sub 
% structure. 
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zercs bb,idi 


r.zc = ,r..ici,c:; zeros(bb,dd. ; 

zeros:dd,bb) h2:c2,o; !; 

% 

hce= ;hici,ili zerosibb.oc) hiicl,cl) 
zeros'ad,aa) h2io2,i2, zeros'dd,bb) 

% 

% We can now perform the synthesis: 

zr = (k + j»comega*l<*exp(-q*Comega) ) * eye(22) 

hccr = M' • hCC * M; 

heescar = hee - hec ♦ M * inv(inv(zr) + hccr i 

HH(counc,=r.eescart8,8) ; 

end; 

stime (clocJc, cO) ,flops 
% 

Comega=.l:.1:50; 

Freq=Comega/(2*pi); 

plot(Freq,2 0*log(HH)) , grid 

xlabeK' Frequency Hz ') 

ylabel('FRF at coordinate of interest dB ') 


% FRF_lndsprdam_data7 

% 

% This is the data file for the synthesis program. 

% The fol-owing data will be provided by this file. 

% 

% internal coords sub structure 1 

% i2- internal coords sub structure 2 
% cl- coruiection coords sub structure 1 
% c2- connection coords sub structure 2 
% 

ii=(l 2 3 4 5 6 7 8 9 16 17 18 19 20 21 22 23 24 25 26 27 28 29 301 ; 

Cl=[10 11 12 13 14 15 31 32 33 34 35 36); 

i2=(4 5 6 7 8 9 10 11 12] ; 

C2=[l 2 3 13 14 15 16 17 18 19 20 21] ; 

% 

% the following is the mapping matrix 
% the mapping matrix is not general and is 
% case specific 
% 

M=[eye(12) ; 

000000 - 100000 ; 

0000000 - 10000 ; 
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00 OOOOOO -ICOO; 
000 - 100000000 ; 

0 i 00-10000000,• 
o.■;coo-lOooooo,• 
0:c000000-100; 
0000000000 - 10 ; 

0 0 0 0 D 0 0 0 0 0 0 - 1 ; 
-lC■000:000000; 
0 - 10000000000 ; 
00-100000000 0 ); 

% spring constant and damping coefficient 
25; 
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APPKVDIX G MATLAB CODE FOR EXAMPI.E SEVEN 


% unif ir-.oelscrsss 

% 

^ This program will calculate the eigenvalues, eigenvectors 
% natural frequencies, and stress frequency response function matrix 
% for a three degree of freedom at each node element. 

% The system is modeled with beam elements that are 
% aligned in the same plane but at any angle (2-D). 


% This program works for a beam element modeled with six general 
% coordinates and thus six DOF. 


% 

% the user must enter the following data to meet the beam configuration 

% (E) youngs modulus psi 

% (I) area moment of inertia in"4 

% (WTD) weight density lbf/in*3 

% (A) cross sectional area in"2 

% conductivity [ the node connection mapping 1 

% node coordinates [ cartesian coordinates for each node 1 

% (beam)the beam element of interest 

% (cc) structure connection coordinates 

% (bcoord) beam element coordinates 

% (cd) distance from beam center to outer most fiber 

% call Che data file 
hullscressdaca 


% calculate the number of beam elements 

% calculate the number of nodes. 
b=size(coord); 
nodes=b(l); 

% 

% convert Che coordinates in to the correct units (in.) 
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% dpply the boundary conditions 

% Che user xusc adjust the giocai .racrix to meet the boundary conditions 
% 

% to delete rows 
kg iBC;, : 1 = .: ; 
mg I .51 , ; = ; 

% to delete tc.umns 
kgi :, IBC. ) = ; 

% 

% 'all the function and calculate eigenvectors and the eigenvalues 
% which are the mode shapes and natural frequency in (rad/seci‘2 

, .ambda,phi] =fgTr.odes (kg, mg) ; 

% this converts the eigenvalues to nac frequency in (rad/sec) and 
% and hertz(1/sec 1 

omega - sqrc(lambda); 
freq = omega/',2*pi) ; 

% 

% calculate the FTIF over the freq. of interest 
* 

counc=0; 

for Comega=2:1:1600 
count=counc+i; 

Z=kg-Comega"2 *mg; 

H=inv(Z); 

% save the coordinate of interest to plot 
HH (count) =H (i, ; 

% 

% this portion calculates the stress FRF in a given element 


HCOL=size(cc); 

NHCOL=HCOL(2); 

MEQ=[0 l(beam)/2 1 0 0 0],- 
[Trmatrixl=fcrans(c,s,beam); 

[kell -fkelemencS(1(beam),WTD(beam),I(beam),E(beam),A(beam),r(beam)); 
for i=l:NHCOL 
HEL-H(:,cc(i)); 

HELR=HEL(bCOOrd,:); 

HLOCAL,=Trmat r ix'HELR ; 

NODEF=kel*HLOCAL; 

stress(1,i)=cd/l(beam)* meq*nodef; 

Stress(counC)=stress(l,3); 

% this is the end of the stress FRF calculation 

Comega=2:1:16 0 0; 

Freq=Comega/(2*pi); 
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pl^c!Freq,20*log(HH)■,grid,pause 
ploc(Freq,20*log(Stress' ,grid 



% r.ullsCressdata 

% 

% 

% T'nis is Che data for che finite element program. 

% che data will be in the form of 
% - (E) youngs modulus psi 
% - (I) area moment of inertia in*4 
% - (WTD) weight density Ibf/in* - 
% - (A) cross sectional area in'. 

% - conductivity [ the node conr. = rcion mapping ) 

% - node coordinates [ cartesiar joords for each node 
% c.ne main program will convert c 

% - (beam) che beam element of interest 
% - (bcoord) che beam element coordinates 
% - (cc) structure connection coordinates 
% - (cd) distance from beam center to outer most fiber 
% 

E=[30 30 30 30 30 ]*le6; 

I=[.1666 .1666 .1666 .1666 .1666]*le-3;%l/l2bh*3 b=2, 
A= t . 2 .2 .2 .2 .2] ; 

WTD=[.2832 .2832 .2832 .2832 .28321; 
beam=4; 

CC=[4 5 6 10 11 12] ; 
bcoord=[l 2 3 10 11 12); 

% 

COn=[1,2; 

2,3; 



coord=[10,0; 

20 , 10 ; 

10 , 20 ; 

0,101; % this is in (in.) 

% 

BC=[] ; 


hullstressdatal 


% This is che data for the finite element program. 
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% This is Che outer structure. 

% 

% the data will be in the form of 

% - :e. youngs modulus psi 

% - 'I, area moment of inertia in"4 

% - (wrDi weight density lbf/in"3 

% - A; cross sectional area in*2 

% - conductivity [ the node connection mapping i 

% - node coordinates [ cartesian coords for each node in ft. ] 

% the main pregram will convert to in.. 

% 

E=130 30 30 301*le6; 

x=L.i666 .1666 .1666 .1666) * le - 3 ; %l/12cr." 3 b-2 , h=0.1 
A=i.2 .2 .2 .2]; 

WTr)=i.2832 .2832 .2832 .2832.; 

% 

OOn=[1,2; 


'’ccrd=[10,0; 

20 ,- 0 ; 

10 , 20 ; 

0,10]; % this is in (in.) 

% 

BC=11 ; 


% hullstressdaca2 

% 

% 

% This is the data for the finite element program. 

% This is Che inner structure. 

% 

% Che data will be in the form of 

% - (E) youngs modulus psi 

% - (I) area moment of inertia in'4 

% - (vmi) weight density lbf/in"3 

% - (A) cross sectional area in*2 

% - conductivity [ the node connection mapping ) 

% - node coordinates [ cartesian coords for each node in ft. 
% the main program will convert to in.. 

% 

E=[30)*le6; 

I=[.1666 ]♦le-3;%i/12bh*3 b=2, h=0.1 
A=[.21 ; 

WTD=[.2832]; 


, 2 ] ; 
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% this 


% program FRF_scress8 

% '-.-.is program will calculate stress FRF by the indirect 
% ''cupliog method 


%lcad exSa.mat % kl, ml are stored here 
%load exSb.mat % k2, m2 are stored here 


stressdataS 

% We need to partion the H matrix of the structure to 
% be modified in the following way. 


% - 



% [hee) = [ ii ic 

% ci cc 1 

% [hec] = [ ic 

% cc 1 

% (hcc) = [ cc ] 

% 

% define the variables 

% internal moment equation in vector form 
MEQ= [0 11/2 1 0 0 0); 

% number of columns in the frf matrix 
nhcol=nodes'*3 ; 

% calling functions 

[kel]=frfkelement6(ll.WTD,I,E,A,rr); 

[TrmatrixJ =frftrans(ccc,ss) 

% set Che clock and flops 

t0=clock; 

flops (0) ; 

counc=0; 

for Comega=2:1:1600 
% 

% Form Frequency Response Models for Each Substructure 

% ----- 
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"ncrdinace Particioning 


z=k2 -"omega * 2 *m2; 


% 

% 



% 


% caioulate the stress of desired beam 
for i=i:nhcol 

heir=helibcoord,: i ; 
hlocal=Trma t rix * he1r; 
nodalF=kel*hlocal; 

S(1,i)-cd/I*MEQ*nodalF; 
end; 

hse=[Sil, ic) S(i,cc) i ; 

% now we will synthesize the stresses 
hsestar=hse - hsc * inv inv(z) * hoc) * hce; 
HS(countl = hsestar(1,9); % this is for coord 6 

end 

etime(clock, to) ,flops 
Comega=2:l:1600,- 
Freq=Comega/(2*pi); 
plot(Freq,20*log(HS)) ,grid 
xlabel(’ Frequency Hz ') 

ylabel('Stress FRF at coordinate of interest dB ’) 


% stressdataS 

% This is the data file for the stress synthesis program. 
% The following data will be provided by this file. 

% 

% ic- internal coords of synthesized structure 
% cc- connection coords of synthesized structure 


% 



7 8 91 ; 

10 11 121 ; 


bcoord=[l 2 3 10 11 12) ; 

% 

B=1301'lee; 

I=[.1666 1*le-3;%l/12bh-3 b=2, h=0.1 


.21 ; 


WTD=[.2832) ; 


cd=.05i 


nodes=-4; 
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APPENDIX H MATLAB CODE EOR EXAMPLE EIGHT 





% unifineeli 

% 

% This program will calculate the eigenvalues, eigenvectors 
% natural frequencies, and frequency response function matrix 
% for a three degree of freedom at each node element. 

% The system is modeled with beam elements that are 
% aligned in the same plane but at any angle {2-D). 


% 


% This program works for a beam element modeled with six general 
% coordinates and thus six DOF. 


% 

% the user must enter the following data to meet the beam configuration 
% (E) youngs modulus psi 
% (I) area moment of inertia in"4 
% (WTD) weight density lbf/in*3 
% (A) cross sectional area in"2 
% conductivity [ the node connection mapping 1 
% node coordinates [ cartesian coordinates for each node 1 
% (bb) structural proportional damping constant 

% call the data file 
fxbeam_data9 

% start the program clock and flops to determine program running 
% time and floating point calculations 
tO=clock,- 
flops(0); 

% calculate the number of beam elements 
a=size^con); 
nume1=a(1); 

% 

% calculate the number of beam elements proportionally damped 
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% call crig 
ic,sl=ftrig( 
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equal damped matrix 


% calculate the beam element lengths and beam angles 
% for the damped beams 
for i=l:numel_damp 

ID=dcon(1,2); 

1 , i) =sqrt ( (coord(ID, 1) - coord (IC, 1) ) "2 + (coord(ID, 2) 

DX(i)=coord(ID, 1) -coorddC, 1^ ; 

DY ( i) =coord (ID, 2) - coorddC, 2) ; 
if DX(i)>=0 & DY(i)>=0; 

t (1, i)=acos(DX(i) /Id, i) ) ; 
elseif DX(i)<0 & DY(i)>=0; 

t (1, i)=acos (DY(i) / I (1, i) ) ♦pi/2,- 
elseif DX(i)<0 & DY(i)<=0; 

t(1,i)=acos(abs(DX(i))/l(l,i))+pi; 


t (1, i)=acos (abs (DY(i))/l(l,i>) + (3*pi/2) ; 


end,- 

% 

% call trig function 

[c, s] =f trig (t, numel_danTp) ; 

% 

% calculate radius of gyration 
for i=l:numel_damp 
r(1,i)=sqrc(I(i)/A(i)); 


for u=l:numel_damp 

[IteU =f element 6 (1 (u) , wtd(u) ,I (u) ,E(u) ,a(u) ,r(u) ,c(u) ,s(u) ) ; 

v=dcon(u,1); 
w=dcon(u,2); 


kgd(3*v-2:3*v,3*v-2:3»v) 
kg<i(3*v-2 ; 3*v, 3*w-2 : 3*w) 
ltgd(3»w-2 : 3*w, 3*v-2 :3*v) 
Icgd ( 3*w - 2 : 3*w, 3*w-2 : 3*w) 
% 

end 


= Icgd (3*v-2 : 3*v, 3*v-2 : 3*v) 
- )cgd(3*v-2:3*v,3*w-2:3*w) 
= Itgd (3*w-2 : 3*w, 3»v-2 : 3*v) 
= kgd(3*w-2:3*w,3*w-2:3*w) 


j*bb*kel 
j »bb*kel 
j•bb*kel 
j*bb*kel 
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% now apply the boundary conditions 

% Che user must adjust the global matrix to meet the boundary condici 

% t ; delete rows 
kg' [BCl , : : = il ; 
kgd '. 3Cj , : ' = :i ; 

% CO delete columns 
kg', (BCi ) = [1 ; 
kgd l :, [BCj = [1 ; 
mg(: , [BCI ) = il ; 

% 

% 'tall the function and calculate eigenvectors and the eigenvalues 
% which are the mode shapes and undamped natural frequency in 
(rad/sec)"2 
% 

. lambda, phi J “fgniodes (kg, mg) ,- 

% this now converts Che eigenvalues to nac frequency in (rad/sec) and 
% hertz(1/sec) 

omega - sqrc(lambda); 
freq = omega/(2*pi); 

% 

% conscuct Che frequency response plot over Che frequencies 
count=0; 

for Comega=.1:.5:500 
counc=count+1; 

Z=kgd-Comega'2 »mg; 

H=inv(Z) ,- 

% this determines the coordinate of inCrest to plot 

HH(count)=H(2,6); 

end; 

% end Che program clock and flops 
ecime(clock,to),flops 
Comega=.i:.5:500; 

Freq=Comega/(2*pi); 

plot(Freq,20*log(HH)),grid 

xlabel(' Frequency Hz ') 

ylabel('FRF at coordinate of interest dB ') 




% unifineelimode 

% This program will calculate Che eigenvalues, eigenvectors 
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% coordinates and 


% Che user must enter the following data to meet the beam configuration 
% (S) youngs modulus psi 
% (I) area moment of inertia in*4 
% (WTD) weight density Ibf/in'l 
% (A) cross sectional area in*2 
% conductivity [ the node connection mapping ) 

% node coordinates [ cartesian coordinates for each node 1 
% (bb) structural proportional damping constant 
% 

% call the data file 
fxbeam_data9 

% start Che program clock and flops to determine program running 
% time and floating point calculations 
cO=clock; 
flops(0); 

% 

% calculate the number of beam elements 
numel=a(i); 


% calculate the number of beam elements proportionally damped 
aa=size(dcon); 
numel_damp-aa(1); 

% calculate the number of nodes. 
b=size(coord); 
nodes=b{l) 

% 

% convert the coordinates in to the correct units (in.) 
coord=coord*i2; 

% calculate the beam element lengths and beam angles in radians 
for i=l:numel 


ID=con(i,2); 
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i i i, i) =sqrc ( (coord113, P - coord - IC, 1. ; "2-' icoord;ID, 2 ■ -coord' :c ,2 . 
DX(i)=coord(ID,i)-coord 12,1'; 

DYii;=coord(ID,2)-coord;iC,2); 
if DXIi'>=0 i DY(i)>=0; 

c ( i, i =acos (DX ( i) /111, i) ) ; 
elseif DX(il<0 & DY(i)>=0; 

c ., i ' =accs 'DY i) /I - , i . ■^pi/2 ; 
elseif DX ii <:■ i DY (i <=' ; 

c(i,i)=acos(abs(DXii !/Ii1,i))*pi; 
else 

t (1, i ; =acos (abs (DY( i - /i (1, ii - * i3*pi/2) ,- 
end; 
end; 


% call trig function 
[c,si=ftrig(c,numel); 

% calculate radius of gyration 
for i=l:numel 
r (1, - =sqrc (I (i) /A ! i1 ; 

% 

% create Che global matrix which is all zeroes. 


kg=[zeros(nodes*3,nodes*3)); 
mg=[zeros(nodes*3,nodes*3)1; 

% 

% assemble the elemental matricies to the global matrix. 


kg(3*v 
kg (3*’ 
kg(3*w 
kg(3*w 

mg(3*v 
mg(3*v 
mg(3 *w 
mg(3 *w 


3*v) = kg(3*v 
3*w) = kg(3*v 
3*v) = kg(3*w 
3*w) = kg(3*w 

3*v) = mg(3*v 
3*w) = mg(3*v 
3*v) = mg(3*w 
3»w) = mg(3*w 


3*v) » keld 
3*w) keld 
3*v) 1- kel(4 
3*w) kel (4 

3*v) meld:3 
3*w) meld:3 
3*v) + mel(4:6 
3*w) t mel(4:6 


apply structural prop, damping to the k matrix and set global 
k matrix to equal damped matrix 
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% calculate the beam element lengths and beam angles in radians 
% for the damped beams 
for i=i:numel_damp 
IC=dcon(i,1); 

ID=dchn!i,2); 

1(1,i)=sqrt((coord{ID,i)-coord IC,1))*2*(coord(ID,21-coord(IC,2))*2); 
DX( ii=.-:oord(:D, I) -coorddC, 1) ; 

DY .i ID,2; -coorddC,2) ; 

if DX,i)>=0 & DY(i)>=0; 

- d, il=acos(DX(i) /I (1, i) ) ; 
elseif DX(i)<0 & DY(i)>=0,- 

t U, i)=acos(DY(i) /I (1, i) j -pi/2; 
elseif DX(i)<0 & DY(i)<=0; 

c (1, i)=acos (abs (DX ' i; ) /I (1, i) ) -pi ,- 
else 

t(1,i)=acos(abs(DY(i))/I(l,i))- (3*pi/2) ; 

end; 

% 

% call trig function 
[c,si=fcrig(c,numel_damp); 

% 

% calculate radius of gyration 
for i=l:numel_damp 
r(1,i)=sqrt(I(i)/A(i)); 
end; 

% 


)cgd=kg; 

for u=l:numel_damp 
[kel]=felements(1(u),WTD( 
% 

v=dcon(u,1); 
w=dcon(u,2); 

% 

kgd(3*v-2:3*v,3*v-2:3*v) 
kgd(3*v-2:3*v,3*w-2:3*w) 
kgd(3*w-2:3*w,3*v-2:3*v) 
kgd(3*w-2:3*w,3*w-2:3*w) 


I ,I (u) ,E(u) ,A(u) ,r(u) ,c(u) ,s{u) ) ; 


kgd(3*v-2:3*v,3*v-2:3*v) + j*bb*kel(1:3,1:3) 
kgd(3*v-2:3*v,3*w-2:3*w) - j*bb*kel(1:3,4:6) 
kgd(3*w-2:3*w,3*v-2:3*v) - j*bb*kel(4:6,1:3) 
kgd(3*w-2:3»w,3*w-2:3*w) - j*bb*kel(4:6,4:6) 


% apply the boundary conditions 

% the user must adjust the global matrix to meet the boundary conditions 

% to delete rows 
kg( [BCl , :) = (1 ; 
kgddBCl,:) - [1 ; 
mg( [BC] , : ) = [] ; 

% to delete columns 
kg(:, [BCl ) = [) ; 
kgd(:, [BCl ) = [] ; 
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% call Che funccion and calculace eigenveccors and che eigenvalues 
% which are che mode shapes and undamped natural frequency in 
(rad/sec)‘2 
% 

:lambda,phi]=fgmcdes!kg,mgl; 

% this new converts che eigenvalues to nat frequency in rad/seo 
% hertz(1/sec) 

% 

omega = sqrc(lambda); 
freq = omega/(2*pi); 

% conscucc the frequency response plot over che frecfuencies 

% of interest 

counc=0,- 

fs=size(freq); 

nm=fs , I; ; 

% extract the mode shapes for the coordinates of interest 
phired=phiici,:); 
for Comega=.1:.5:500 
counc=counc+i; 

% generate che diagonal frequency matrix 
for i=i:nm 

nomega{i)=1/(omega(i)"2-Comega"2); 

% generate che FRF matrix 
H=phired»diag(nomega)'phired'; 

% this determines che coordinate of increst to plot 
HH(counc)-H(i,3); % this refers to (2,6) 

end; 

% end the program clock and flops 
etime(clock, to) ,flops 
Comega=.i:.5:500; 

Freq=Comega/(2*pi); 

plot(Freq,20»log(HH)) ,grid 

xlabel(’ Frequency Hz ') 

ylabel('FRF at coordinate of interest dB ') 


% fxbeam_daca9 


% This is che data for the finite element program and che finite 
% element program using modal representation. 

% 

% the data will be in che form of 
% - (E) youngs modulus psi 

% - (I) area moment of inertia in*4 
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% - iWTD) weight density lbf/in-3 
% - !A) cross sectional area in -2 
% - conductivity [ the node connection mapping ] 

% - node coordinates i cartesian coords for each node in ft. 1 
% the main program will convert to in.. 

% - ibbj structural proportional damping constant 
% - ri' 'ccrdinates of interest 
E=;3o ii ; 7; ; *ie6; 

I=[.02083 .02083 .02083 .02083] ;%l/12bh-3 b=2, h=.5 
A= : 1 1 i 1 i ; 

WTD=[.2832 .2832 .2832 .2832]; 
bb=0 ; 

% coordinates of interest (any internal and all connection) 
ci=[2 4 5 61 ; 


coord=[0, 

8 , 


,0 

, 0 


12 , 0 ; 

16,0] ; 


BC=[1 2 3 13 14 15]; 


% right_struc9 

% 

% This is the data for the finite element program. 

% 

% the data will be in the form of 

% - (E) youngs modulus psi 

% - (I) area moment of inertia in'4 

% - (WTD) weight density lbf/in'3 

% - (A) cross sectional area in'2 

% - conductivity [ the node connection mapping ] 

% - node coordinates ( cartesian coords for each node in ft. ] 
% the main program will convert to in.. 

% • (bb) structural proportional damping constant 
% 

E=[30 30]*le6; 

I=t .02083 .02083];%l/12bh'3 b=2, h=.5 
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bb=0; 

% 

2,31 ; 

% 

dcon= L1,2; 
2,3; ; 


coord= 10,0; 



BC=17 8 9]; 




% This is the data for the finite element program. 

% the data will be in the form of 

% - (E) youngs modulus psi 

% - (I) area moment of inertia in"4 

% - (WTD) weight density lbf/in*3 

% - (A) cross sectional area in'2 

% - conductivity [ the node connection mapping ] 

% - node coordinates ( cartesian coords for each node 
% the main program will convert t 

% - (bb) structural proportional damping constant 
E=[30 30)*le6; 

I=[ .02083 .02083];%i/i2bh*3 b=2, h=.5 
WTD=[.2832 .2832] ; 



2,3] ; 


dcon=[l,2; 
2,3] ; 

coord=[0,0; 

4,0; 

8 , 0 ]; 

% 

BC=[1 2 3] ; 


Clear 

cig 

% FRF_Synth9 
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% nil= redefined internal coords of sub structure 1 
% ncl= redefined connection coords of sub structure 1 
% ni2= redefined internal coords of sub structure 2 
% nc2= redefined connection coords of sub structure 2 
% 

FRF_Synth_data9 



% extract the rows of phi relating to the coordinates of interest 
phiredl=phil(cii,:); 
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phired2=phi2(ci2, :) ; 
for coniega=. i:. 5:500 
couiit=count 1 ; 

% 

% Form Frequency Response Models for Each Subscruccure 


for i=::5 

r.onega- i = - / omega 1' i i * 2 - Comega" 2) ; 
end ; 

hl=phiredl*diag(nomegal)‘phiredl'; 
for i=i:6 

nomega2(i)=1/iomega2(i)"2-Comega‘2); 
h 2 =phlred 2 *diag(nomega 2 )*phired2'; 



cc-c(2) 
dd=d (2) 


Remember, 




calculate the following: 




So we need to assemble [hie] and [hcc] using the 
the coordinate sets we just defined. 

These matrices contain the FRF data for both substructures 
prior to coupling, i.e the pre-synthesis FRF data. 


Coordinate Partitioning 


Build up uncoupled FRF matrix and sub-partitions: 


% 

hie = [hl(nil,ncl) zeros(aa,dd)]; 


hcc = [hl(ncl,ncl) zeros(bb,dd); 

zeros(dd,bb) h2(nc2,nc2)]; 

% 

% We can now perform the synthesis: 
hccr=M' * hcc * M; 

hicstar = hie - hie * M * inv( heer ) * M' * hcc; 


% look at the coordinate of interest 
HH(count)=hicstar(1,3); 
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FHF_Synth_dai:a9 

% 

c ^^rlock; 
flops-O ; 
counc=0; 

% exr.rac- the rows of phi relating to the coordinates of interest 

phiredi=p;".i - icii, : ) ; 

phired2=phi2(ci2,:); 

for Comega=.l:.5:500 

counc=count+l; 

% 

% Form Frequency Response Models for Each Substructure 


for i=l:6 

nomegal(i)=1/(omegal(i)*2-Comega‘2); 


hi=phiredl*diag(nomegal)'phiredl'; 
for i=L:6 

nomega2(i)= 1 /(omega2(i)•2-Comega'2); 
h2=phired2*diag(nomega2)*phired2'; 


b=size(cl) 



dd=d(2) 


% Remember, we are trying to calculate the following-. 

% hccr = M’ ♦ hcc * M 
% 

% So we need to assemble [hcc] using the 
% the coordinate sets we just defined. 

% These matrices contain the FRF data for both substructures 
% prior to coupling, i.e the pre-synthesis FRF data. 


% Coordinate Partitioning 

% -- 

% Build up uncoupled FRF matrix and sub-partitions: 
% 

hcc = (hi(ncl,ncl) zeros(bb, dd); 

zeros(dd,bb) h2 (nc2 ,nc2) 1 ,- 

% 

% We can now perform the synthesis: 
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% determine the nat freq of the synthesized structur 
dhccr=det(hccr); 
dh!count)=dhcc r; 

etime(clock,tO).flops 
Comega=.1:1:500; 

Freq-Conega/!2*pi;; 

axis( iO 80 -le-12 le-12i) 

plot(Freq.dh).grid 

xlabel’' Frequency Hz ') 

ylabel('Determinent of Hoc reduced') 
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APPENDIX I GENERAL MATLAB FUNCTIONS 



function [kel,mel] = felements(1,WTD,I,E,A,r,c,s) 

% This function is generating the elemental mass and 
% stiffness matrix. Input is element length, weight density, 

% area moment of inertia, young modulus, radius of gyration 
% and the angel of the element with respect to horizontal x 
% axis. Counter clock is positive angle and clock is negative % angle. 



mel(2,1) 
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element 








kel (. 


,<ei 

kel(3,4i 
kel(1,1) 
kel( 1,2) 
kel(4,3) 
kel (1,5) 
kel(2,5) 
kel(3,5) 
kel(4,5) 
kel(2,2) 
kel(3,5) 
kel(1,6) 
kel(2,6) 
kel(3,6) 
kel(4,6) 
kel(5,6) 
kel(3,3) 


% now apply cine stiffness constant to the elemental matrix 


kel=kel*(E*I/(1*3)); 


function[lambda,phi]=fgmodes(kg,mg) 

% This function calculates the natural frequencies and 
% the mode shapes. Theses are Che eigenvalues and 
% eigenvectors 
% 

a=length(mg); 

[V, D] =eig(mg\kg); 

(omga,index)-sort(diag(D)); 
lambda=zeros(a,a); 
for i=l:a; 

lambda(i,i)=omga(i); 


phltemp ( : , i) =v (:, index (i) ) ,■ 
end,■ 

lambda=diag(lambda); 

[phi,orth]=fgmassnorm(phicemp,mg); 
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kdsprg = zeros(6); 
kdsprg (1,1) = 1< + 
lcdsprg(l,4) = -k - 
kdsprg(2,2) =k* 
kdsprg (2,5) = -k - 
kdsprg (3 , 3) = k 
kdsprg(3,6) = -k - 
kdsprg (4,1) = -k - 
kdsprg(4,4) = k + 
kdsprg(5,2) =-k- 
kdsprg(5,5) =k+ 
kdsprg(6,3) =-k- 
kdsprg(6,6) = k 


j*Cotnega*B*k; 

j»Comega*B*k 

j*Comega*B»k; 

j*Comega*B*k 

j»Comega*B*k; 

j*Comega»B»k 

j*Comega*B»k 

j*Coinega*B*k,- 

j»Comega*B*k 

j*Comega*B*k; 

j »Comega*B*k 
j•Comega»B*k; 


funct ion [kdsprg] =f sprngdainpC (k, Comega, q) ; 


% This function will generate a spring and damper system k 
% matrix of size [6 x 6] which correlates with the 3 dof of 
% the beam element kx, ky, k(theta), cx, cy, c(theta). 

% Remember that for a spring of stiffness k and a damper with 
% damping (Co e"-qO), the matrix looks like 
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Trmacrix(5,5 > =ccc; 
Trmat;rix(6,6) =i ; 


179 






LIST OF REFF.RENCFS 


1. Gordts. J.H..Bielawa. R.L.. Flannelly, W.G.. "A General Theory for Frequency 
Domciin Structural Synthesis,” Journal of Sound and Vibration, Vol. 150 No. 1. pp. 
139-158, 1989. 

2. Gordis. J. FI., Flannelly, W.G. Analysis of Stress Due to Fastener Tolerance in 
Assembled Components. Proceedings of the 34th AIAA/ASMII/ASCE/AHS/ACS 
Structures, Structural Dynamics, and Materials Conference, La Jolla, CA., 1993. 

3. Gordis, J. FI. Structural Synthesis in the Frequency Domain: A General 
Formulation. To appear. Proceedings of the 12th International Modal Analysis 
Conference, Honolulu, HI., 1993. 

4. James, M.L., Smith, G.M., Wolford, J.C., Whaley. P.W., Vibration of Mechanical 
and Structural Systems. Harper and Row, Publishers, Inc., New York, NY. pp. 
482-572. 1989. 


180 





INITIAL DISTRIBUTION LIST 


No. Copies 


1. Defense Technical Information Center 2 

Cameron Station 

Alexandria, Virginia 22304-6145 

2. Library, Code 52 2 

Naval Postgraduate School 

Monterey, California 93943-5002 

3. Professor!. H. Gordis, Code ME/Go 3 

Department of Mechanical Engineering 

Nava] Postgraduate School 
Monterey, California 93943 

4. Professor A. Healey, Code ME/He 1 

Department of Mechanical Engineering 

Naval Postgraduate School 
Monterey, California 93943 

5. Department Chairman, Code ME 1 

Department of Mechanical Engineering 

Naval Postgraduate School 
Monterey, California 93943 

6. Naval Engineering Curricular Office (Code 34) 1 

Naval Postgraduate School 

Monterey, California 93943 

7. LT. Ronald E. Cook 1 

1090 Derby Street 

Casper, Wyoming 82609 


181 








OiJOi-E Y KNOX LIBRARY 
lAYAL POSTGRADUATE SCHOOl 
MONVeREY CA 93943-5101 









3 2768 00019545 7 






